ABCD Media GFA Suicide Analysis 08.31.2018

myGAMM4Bin <- function(dv,iv,cv,nv,dat)
{
  
  indv <- paste(iv, collapse=" + ")
  cova <- paste(cv, collapse=" + ")
  if(length(nv)>1){
    nstv <- paste("~","(","1","|",nv[1],"/",nv[2],")",sep="")
  } else {
    nstv <- paste("~","(","1","|",nv[1],")",sep="")
  }
  datnames <- names(dat)
  
  if(iv %in% datnames) {
    form1 <- paste(dv," ~ ",indv," + ",cova,sep="")
  } else { form1 <- paste(dv," ~ ",cova,sep="")}
  
  # print(form1)
  # print(nstv)
  
  mygam <- gamm4(as.formula(form1), family=binomial, random = as.formula(nstv), data = dat)
  
  return(mygam)
}

mydir <- paste0("/Users/mpaulus/Dropbox (Personal)/Private/RDataAnalysis/ABCD_Data/Media/")
# File components:
myfile <- c("ABCD_")
datatext <-("SMA_CBCL_COG_ENV_")
GFAtext <- c("R_GFA_")
dateext <- c("09.02.2018")

myall <- paste(mydir,myfile,datatext,GFAtext,dateext,".RData",sep="")
load(myall)
mynames <- names(abcdnegreinf)

# Load different variable sets:

myall <- paste(mydir,"abcd_activars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_cbclvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_cogvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_friendvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_medvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_sulcvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_thickvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_volvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_screenvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_socialsummaryvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_physvars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_covars_08.23.2018",".RData",sep="")
load(myall)
myall <- paste(mydir,"abcd_suicidvars_08.23.2018",".RData",sep="")
load(myall)

# Assigns a temporary data set for computational purposes:
currdata <- abcdnegreinf

# Rename kids suicide variables:

kidsuicidevars <- suicidvars[intersect(grep("_t",suicidvars),grep("ksads",suicidvars))]

kidsuicidelabels <- c("selfinjurious_now","selfinjurious_past","wishdead_now","wishdead_past",
                      "SI_present","SI_past","suicideattempt_now","suicideattempt_past",
                      "selfinjurytodie_now","selfinjurythought_now","suicidethoughtmethod_now",
                      "suicideintent_now","suicideplan_now","suicideprep_now","suicideabort_now",
                      "suicidemethod_now","suicideattempt_now","selfinjurytodie_past",
                      "selfinjurythought_past","suicidethoughtmethod_past","suicideintent_past",
                      "suicideplan_past","suicideprep_past","suicideabort_past",
                      "suicidenumber_past","suicidemethod_past","suicideattempt_past")

# Rename the variables:
setnames(currdata, old=c(kidsuicidevars), new=c(kidsuicidelabels))

# Rename parents suicide variables:
parsuicidevars <- suicidvars[intersect(grep("_p",suicidvars),grep("ksads",suicidvars))]
parsuicidevars <- parsuicidevars[-c(1:4)]

parentsuicidelabels <- c("p_suicideAttempt_now","p_suicideAttempt_past","p_selfinjury_now","p_selfinjury_past","p_SImethod_now","p_suicideintent_now","p_suicidePlan_now","p_suicidePrep_now","p_suicudeInterrupt_now","p_suicidemethod_now","p_suicideTought_now","p_selfinjuryDeath_past","p_selfinjuryThought_past","p_SImethod_past","p_suicideintent_past","p_suicidePlan_past","p_suicidePrep_past","p_suicideInterrupt_past","p_numSA_past","p_Samethod_past","p_expectToDie_past")

homicidallabels <- c("p_homicidaldeas_now","p_homicidaldeas_poast","p_homicidalPlan_now","p_homicidalPlan_past")

# Rename the variables:
setnames(currdata, old=c(parsuicidevars), new=c(parentsuicidelabels))

kidsuicidetotal <- rowSums(currdata[,kidsuicidelabels])
parsuicidetotal <- rowSums(currdata[,parentsuicidelabels])

currdata$kidsitotal <- c(kidsuicidetotal)
currdata$parsitotal <- c(parsuicidetotal)

# Create a binary variable
currdata$KidsSIyes <- ifelse(currdata$kidsitotal>0,1,0)
currdata$ParSIyes <- ifelse(currdata$parsitotal>0,1,0)
currdata$KidsSIyes <- as.factor(currdata$KidsSIyes)
currdata$ParSIyes <- as.factor(currdata$ParSIyes)

# Describe and characterize the total variables
# individual kids items:
factordata <- data.frame(lapply(currdata[,kidsuicidelabels], factor))
describe(factordata[,kidsuicidelabels])

## factordata[, kidsuicidelabels] 
## 
##  27  Variables      4524  Observations
## ---------------------------------------------------------------------------
## selfinjurious_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4341   130
## Proportion 0.971 0.029
## ---------------------------------------------------------------------------
## selfinjurious_past 
##        n  missing distinct 
##     4471       53        2 
##                     
## Value         0    1
## Frequency  4249  222
## Proportion 0.95 0.05
## ---------------------------------------------------------------------------
## wishdead_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4415    56
## Proportion 0.987 0.013
## ---------------------------------------------------------------------------
## wishdead_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4170   301
## Proportion 0.933 0.067
## ---------------------------------------------------------------------------
## SI_present 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4400    71
## Proportion 0.984 0.016
## ---------------------------------------------------------------------------
## SI_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4315   156
## Proportion 0.965 0.035
## ---------------------------------------------------------------------------
## suicideattempt_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4464     7
## Proportion 0.998 0.002
## ---------------------------------------------------------------------------
## suicideattempt_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4443    28
## Proportion 0.994 0.006
## ---------------------------------------------------------------------------
## selfinjurytodie_now 
##        n  missing distinct    value 
##     4471       53        1        0 
##                
## Value         0
## Frequency  4471
## Proportion    1
## ---------------------------------------------------------------------------
## selfinjurythought_now 
##        n  missing distinct 
##     4471       53        2 
##                     
## Value         0    1
## Frequency  4470    1
## Proportion    1    0
## ---------------------------------------------------------------------------
## suicidethoughtmethod_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4447    24
## Proportion 0.995 0.005
## ---------------------------------------------------------------------------
## suicideintent_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4459    12
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## suicideplan_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4459    12
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## suicideprep_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4462     9
## Proportion 0.998 0.002
## ---------------------------------------------------------------------------
## suicideabort_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4458    13
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## suicidemethod_now 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4451    20
## Proportion 0.996 0.004
## ---------------------------------------------------------------------------
## suicideattempt_now.1 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4464     7
## Proportion 0.998 0.002
## ---------------------------------------------------------------------------
## selfinjurytodie_past 
##        n  missing distinct    value 
##     4471       53        1        0 
##                
## Value         0
## Frequency  4471
## Proportion    1
## ---------------------------------------------------------------------------
## selfinjurythought_past 
##        n  missing distinct    value 
##     4471       53        1        0 
##                
## Value         0
## Frequency  4471
## Proportion    1
## ---------------------------------------------------------------------------
## suicidethoughtmethod_past 
##        n  missing distinct    value 
##     4471       53        1        0 
##                
## Value         0
## Frequency  4471
## Proportion    1
## ---------------------------------------------------------------------------
## suicideintent_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4453    18
## Proportion 0.996 0.004
## ---------------------------------------------------------------------------
## suicideplan_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4459    12
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## suicideprep_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4456    15
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## suicideabort_past 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4454    17
## Proportion 0.996 0.004
## ---------------------------------------------------------------------------
## suicidenumber_past 
##        n  missing distinct 
##     4471       53        2 
##                     
## Value         0    1
## Frequency  4427   44
## Proportion 0.99 0.01
## ---------------------------------------------------------------------------
## suicidemethod_past 
##        n  missing distinct    value 
##     4471       53        1        0 
##                
## Value         0
## Frequency  4471
## Proportion    1
## ---------------------------------------------------------------------------
## suicideattempt_past.1 
##        n  missing distinct 
##     4471       53        2 
##                       
## Value          0     1
## Frequency   4443    28
## Proportion 0.994 0.006
## ---------------------------------------------------------------------------

# total item:
print(describe(currdata$kidsitotal))

## currdata$kidsitotal 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##     4471       53       14    0.347   0.2691   0.4941        0        0 
##      .25      .50      .75      .90      .95 
##        0        0        0        1        2 
##                                                                       
## Value          0     1     2     3     4     5     6     7     8     9
## Frequency   3877   352   114    39    29    24    18     5     4     2
## Proportion 0.867 0.079 0.025 0.009 0.006 0.005 0.004 0.001 0.001 0.000
##                                   
## Value         10    11    12    14
## Frequency      4     1     1     1
## Proportion 0.001 0.000 0.000 0.000

hist(currdata$kidsitotal,main=paste("Histogram of Total Youth Suicid Items"))

# individual parent items:
factordata <- data.frame(lapply(currdata[,parentsuicidelabels], factor))
describe(factordata[,parentsuicidelabels])

## factordata[, parentsuicidelabels] 
## 
##  21  Variables      4524  Observations
## ---------------------------------------------------------------------------
## p_suicideAttempt_now 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4451     3
## Proportion 0.999 0.001
## ---------------------------------------------------------------------------
## p_suicideAttempt_past 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4442    12
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## p_selfinjury_now 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------
## p_selfinjury_past 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------
## p_SImethod_now 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4445     9
## Proportion 0.998 0.002
## ---------------------------------------------------------------------------
## p_suicideintent_now 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4447     7
## Proportion 0.998 0.002
## ---------------------------------------------------------------------------
## p_suicidePlan_now 
##        n  missing distinct 
##     4454       70        2 
##                     
## Value         0    1
## Frequency  4452    2
## Proportion    1    0
## ---------------------------------------------------------------------------
## p_suicidePrep_now 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4449     5
## Proportion 0.999 0.001
## ---------------------------------------------------------------------------
## p_suicudeInterrupt_now 
##        n  missing distinct 
##     4454       70        2 
##                     
## Value         0    1
## Frequency  4453    1
## Proportion    1    0
## ---------------------------------------------------------------------------
## p_suicidemethod_now 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4450     4
## Proportion 0.999 0.001
## ---------------------------------------------------------------------------
## p_suicideTought_now 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4450     4
## Proportion 0.999 0.001
## ---------------------------------------------------------------------------
## p_selfinjuryDeath_past 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------
## p_selfinjuryThought_past 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------
## p_SImethod_past 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------
## p_suicideintent_past 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4441    13
## Proportion 0.997 0.003
## ---------------------------------------------------------------------------
## p_suicidePlan_past 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4451     3
## Proportion 0.999 0.001
## ---------------------------------------------------------------------------
## p_suicidePrep_past 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4437    17
## Proportion 0.996 0.004
## ---------------------------------------------------------------------------
## p_suicideInterrupt_past 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4447     7
## Proportion 0.998 0.002
## ---------------------------------------------------------------------------
## p_numSA_past 
##        n  missing distinct 
##     4454       70        2 
##                       
## Value          0     1
## Frequency   4435    19
## Proportion 0.996 0.004
## ---------------------------------------------------------------------------
## p_Samethod_past 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------
## p_expectToDie_past 
##        n  missing distinct    value 
##     4454       70        1        0 
##                
## Value         0
## Frequency  4454
## Proportion    1
## ---------------------------------------------------------------------------

# Rename some of the variables:
names(currdata)[names(currdata)=="highest.household.income"] <- "HHInc"
currdata$HHInc <- factor(currdata$HHInc,levels = c("[<50K]","[>=50K & <100K]","[>=100K]",""))
levels(currdata$HHInc) <- c("[<50K]","[50K - 100K]","[>100K]","_miss")
levels(currdata$high.educ) <- c("<= 12 grades","HS Degree","College Degree","Bachelor",
                                    "Higher","_miss")

# Creating Quartiles of the first 10 Robust GFAs:
# Alternatively cut by z-scores of the factors < -.5, -.5 - .5, > .5
myqGFA <- paste0("qGFA",rep(1:8))
for(i in 1:length(myqGFA)){
currdata[,myqGFA[i]] <- as.factor(quantcut(currdata[,paste0("SMA_RGFA",i)]))  
}

mynames <- names(currdata)

# General Sample Characteristics:
# Need to rename the demographic variables:
demovars <- c("age","female","race.ethnicity","high.educ","married",
              "HHInc","anthro_1_height_in","anthroweight1lb","demo_prnt_age_v2")

# Additional Screen variables - youth reported hours a week:
addscreenvars <- mynames[grep("week_y_",mynames)]

listvars <- c(demovars,addscreenvars,"kidsitotal","parsitotal","KidsSIyes","ParSIyes")

tabledata <- currdata[,c(myqGFA,listvars)]

catVars = c("female","race.ethnicity","married","high.educ","HHInc","KidsSIyes","ParSIyes")
# length(myqGFA)
for(i in 1: length(myqGFA)){
mytable1 <- CreateTableOne(vars = listvars,data=tabledata,factorVars=catVars,strata=c(myqGFA[i]))

# Make the Table look nice:
tabAsStringMatrix <- print(mytable1, printToggle = FALSE, noSpaces = TRUE) 
print(kable(tabAsStringMatrix, "html",caption = paste0("Quartiles of Robust GFA",i)) %>% kable_styling(bootstrap_options = c( "hover", "condensed", "responsive"),font_size = 11))

}

Quartiles of Robust GFA1

	[-1.51,-0.636]	(-0.636,-0.245]	(-0.245,0.379]	(0.379,5.24]	p	test
n	1131	1131	1131	1131
age (mean (sd))	120.56 (7.18)	119.83 (7.46)	119.94 (7.29)	119.76 (7.40)	0.037
female = yes (%)	606 (53.6)	570 (50.4)	527 (46.6)	449 (39.7)	<0.001
race.ethnicity (%)					<0.001
White	608 (53.8)	696 (61.5)	702 (62.1)	645 (57.1)
Black	151 (13.4)	84 (7.4)	96 (8.5)	109 (9.6)
Hispanic	254 (22.5)	209 (18.5)	199 (17.6)	227 (20.1)
Asian	28 (2.5)	32 (2.8)	28 (2.5)	15 (1.3)
Other	90 (8.0)	110 (9.7)	105 (9.3)	134 (11.9)
high.educ (%)					<0.001
<= 12 grades	51 (4.5)	32 (2.8)	31 (2.7)	62 (5.5)
HS Degree	302 (26.7)	257 (22.7)	271 (24.0)	336 (29.7)
College Degree	153 (13.5)	118 (10.4)	129 (11.4)	157 (13.9)
Bachelor	335 (29.6)	364 (32.2)	370 (32.7)	319 (28.2)
Higher	289 (25.6)	360 (31.8)	329 (29.1)	254 (22.5)
_miss	1 (0.1)	0 (0.0)	1 (0.1)	3 (0.3)
married = yes (%)	817 (72.2)	848 (75.0)	804 (71.1)	736 (65.1)	<0.001
HHInc (%)					<0.001
[<50K]	265 (23.4)	206 (18.2)	240 (21.2)	332 (29.4)
[50K - 100K]	287 (25.4)	330 (29.2)	319 (28.2)	327 (28.9)
[>100K]	472 (41.7)	521 (46.1)	490 (43.3)	378 (33.4)
_miss	107 (9.5)	74 (6.5)	82 (7.3)	94 (8.3)
anthro_1_height_in (mean (sd))	55.54 (3.09)	55.39 (3.18)	55.41 (3.17)	55.37 (3.20)	0.570
anthroweight1lb (mean (sd))	83.66 (23.02)	80.65 (21.97)	81.39 (22.35)	83.52 (23.86)	0.002
demo_prnt_age_v2 (mean (sd))	40.38 (6.81)	40.91 (6.66)	40.58 (6.54)	39.94 (7.09)	0.007
week_y_tvmovie (mean (sd))	8.52 (6.96)	8.71 (7.26)	8.40 (6.95)	9.79 (7.89)	<0.001
week_y_video (mean (sd))	5.54 (6.94)	5.61 (6.96)	6.43 (7.96)	7.35 (8.45)	<0.001
week_y_games (mean (sd))	6.02 (6.99)	6.49 (7.08)	7.12 (7.74)	8.20 (8.38)	<0.001
week_y_text (mean (sd))	1.62 (3.77)	1.31 (3.05)	1.29 (3.21)	1.68 (3.97)	0.011
week_y_socnet (mean (sd))	0.68 (2.35)	0.58 (1.95)	0.61 (2.33)	0.96 (3.48)	0.001
week_y_chat (mean (sd))	1.18 (2.87)	1.03 (2.51)	1.11 (3.00)	1.30 (3.53)	0.153
kidsitotal (mean (sd))	0.14 (0.58)	0.17 (0.76)	0.27 (0.93)	0.50 (1.34)	<0.001
parsitotal (mean (sd))	0.00 (0.04)	0.01 (0.11)	0.01 (0.16)	0.08 (0.51)	<0.001
KidsSIyes = 1 (%)	93 (8.3)	104 (9.4)	143 (12.8)	254 (22.7)	<0.001
ParSIyes = 1 (%)	2 (0.2)	3 (0.3)	5 (0.4)	39 (3.5)	<0.001

Quartiles of Robust GFA2

	[-3.8,-0.637]	(-0.637,-0.0118]	(-0.0118,0.645]	(0.645,3.19]	p	test
n	1131	1131	1132	1130
age (mean (sd))	121.05 (7.29)	120.25 (7.34)	119.42 (7.31)	119.38 (7.29)	<0.001
female = yes (%)	575 (50.8)	568 (50.2)	531 (46.9)	478 (42.3)	<0.001
race.ethnicity (%)					<0.001
White	801 (70.8)	738 (65.3)	659 (58.2)	453 (40.2)
Black	30 (2.7)	64 (5.7)	98 (8.7)	248 (22.0)
Hispanic	134 (11.8)	190 (16.8)	264 (23.3)	301 (26.7)
Asian	47 (4.2)	28 (2.5)	18 (1.6)	10 (0.9)
Other	119 (10.5)	111 (9.8)	93 (8.2)	116 (10.3)
high.educ (%)					<0.001
<= 12 grades	5 (0.4)	19 (1.7)	43 (3.8)	109 (9.6)
HS Degree	152 (13.4)	235 (20.8)	315 (27.8)	464 (41.1)
College Degree	94 (8.3)	115 (10.2)	174 (15.4)	174 (15.4)
Bachelor	395 (34.9)	412 (36.4)	346 (30.6)	235 (20.8)
Higher	484 (42.8)	349 (30.9)	254 (22.4)	145 (12.8)
_miss	1 (0.1)	1 (0.1)	0 (0.0)	3 (0.3)
married = yes (%)	932 (82.4)	881 (77.9)	769 (67.9)	623 (55.1)	<0.001
HHInc (%)					<0.001
[<50K]	123 (10.9)	194 (17.2)	259 (22.9)	467 (41.3)
[50K - 100K]	300 (26.5)	324 (28.6)	345 (30.5)	294 (26.0)
[>100K]	646 (57.1)	539 (47.7)	430 (38.0)	246 (21.8)
_miss	62 (5.5)	74 (6.5)	98 (8.7)	123 (10.9)
anthro_1_height_in (mean (sd))	55.70 (3.02)	55.50 (3.10)	55.27 (3.04)	55.25 (3.45)	0.001
anthroweight1lb (mean (sd))	79.50 (19.58)	81.07 (21.17)	82.40 (23.05)	86.25 (26.47)	<0.001
demo_prnt_age_v2 (mean (sd))	41.90 (5.97)	41.08 (6.05)	40.12 (7.11)	38.70 (7.47)	<0.001
week_y_tvmovie (mean (sd))	6.76 (5.82)	8.25 (6.57)	9.48 (7.54)	10.94 (8.34)	<0.001
week_y_video (mean (sd))	4.44 (6.24)	5.32 (6.70)	6.89 (8.07)	8.27 (8.71)	<0.001
week_y_games (mean (sd))	5.49 (6.45)	6.13 (6.81)	7.49 (7.76)	8.73 (8.81)	<0.001
week_y_text (mean (sd))	0.97 (2.37)	1.13 (2.61)	1.56 (3.57)	2.23 (4.86)	<0.001
week_y_socnet (mean (sd))	0.39 (1.54)	0.51 (2.07)	0.75 (2.49)	1.17 (3.70)	<0.001
week_y_chat (mean (sd))	0.66 (1.53)	0.94 (2.20)	1.12 (2.64)	1.89 (4.59)	<0.001
kidsitotal (mean (sd))	0.25 (0.93)	0.21 (0.80)	0.26 (0.90)	0.35 (1.15)	0.005
parsitotal (mean (sd))	0.01 (0.16)	0.01 (0.19)	0.02 (0.22)	0.05 (0.44)	0.014
KidsSIyes = 1 (%)	152 (13.6)	117 (10.5)	150 (13.4)	175 (15.7)	0.004
ParSIyes = 1 (%)	10 (0.9)	8 (0.7)	12 (1.1)	19 (1.7)	0.130

Quartiles of Robust GFA3

	[-2,-0.564]	(-0.564,-0.193]	(-0.193,0.287]	(0.287,8.66]	p	test
n	1135	1128	1130	1131
age (mean (sd))	119.02 (7.29)	120.01 (7.27)	120.42 (7.34)	120.65 (7.35)	<0.001
female = yes (%)	619 (54.5)	551 (48.8)	507 (44.9)	475 (42.0)	<0.001
race.ethnicity (%)					<0.001
White	695 (61.3)	713 (63.2)	682 (60.4)	561 (49.6)
Black	85 (7.5)	71 (6.3)	95 (8.4)	189 (16.7)
Hispanic	241 (21.3)	212 (18.8)	188 (16.6)	248 (21.9)
Asian	24 (2.1)	32 (2.8)	27 (2.4)	20 (1.8)
Other	89 (7.8)	100 (8.9)	138 (12.2)	112 (9.9)
high.educ (%)					<0.001
<= 12 grades	67 (5.9)	34 (3.0)	33 (2.9)	42 (3.7)
HS Degree	256 (22.6)	243 (21.5)	276 (24.4)	391 (34.6)
College Degree	134 (11.8)	117 (10.4)	133 (11.8)	173 (15.3)
Bachelor	364 (32.1)	373 (33.1)	349 (30.9)	302 (26.7)
Higher	312 (27.5)	361 (32.0)	338 (29.9)	221 (19.5)
_miss	2 (0.2)	0 (0.0)	1 (0.1)	2 (0.2)
married = yes (%)	827 (72.9)	855 (75.8)	826 (73.1)	697 (61.6)	<0.001
HHInc (%)					<0.001
[<50K]	250 (22.0)	217 (19.2)	224 (19.8)	352 (31.1)
[50K - 100K]	306 (27.0)	311 (27.6)	317 (28.1)	329 (29.1)
[>100K]	474 (41.8)	526 (46.6)	505 (44.7)	356 (31.5)
_miss	105 (9.3)	74 (6.6)	84 (7.4)	94 (8.3)
anthro_1_height_in (mean (sd))	55.11 (3.12)	55.42 (3.04)	55.44 (3.06)	55.76 (3.38)	<0.001
anthroweight1lb (mean (sd))	80.31 (22.18)	80.31 (21.84)	81.45 (21.18)	87.15 (25.28)	<0.001
demo_prnt_age_v2 (mean (sd))	40.50 (6.88)	40.75 (6.31)	40.81 (6.72)	39.76 (7.16)	0.001
week_y_tvmovie (mean (sd))	5.54 (4.93)	7.15 (5.73)	9.26 (6.99)	13.49 (8.46)	<0.001
week_y_video (mean (sd))	2.11 (2.76)	3.27 (4.15)	6.06 (6.35)	13.49 (9.48)	<0.001
week_y_games (mean (sd))	2.65 (2.97)	4.04 (4.21)	7.02 (6.33)	14.13 (9.31)	<0.001
week_y_text (mean (sd))	0.34 (0.84)	0.58 (1.01)	1.05 (1.46)	3.93 (6.13)	<0.001
week_y_socnet (mean (sd))	0.08 (0.40)	0.21 (0.65)	0.39 (1.00)	2.14 (4.75)	<0.001
week_y_chat (mean (sd))	0.25 (0.67)	0.52 (1.01)	0.86 (1.44)	2.98 (5.27)	<0.001
kidsitotal (mean (sd))	0.24 (0.91)	0.22 (0.76)	0.25 (0.89)	0.37 (1.19)	0.001
parsitotal (mean (sd))	0.03 (0.31)	0.02 (0.30)	0.01 (0.11)	0.03 (0.32)	0.091
KidsSIyes = 1 (%)	132 (11.8)	135 (12.1)	149 (13.4)	178 (15.9)	0.018
ParSIyes = 1 (%)	16 (1.4)	9 (0.8)	6 (0.5)	18 (1.6)	0.049

Quartiles of Robust GFA4

	[-3.75,-0.47]	(-0.47,0.111]	(0.111,0.528]	(0.528,7.52]	p	test
n	1132	1130	1131	1131
age (mean (sd))	119.92 (7.36)	119.74 (7.35)	120.10 (7.24)	120.34 (7.40)	0.240
female = yes (%)	311 (27.5)	447 (39.6)	617 (54.6)	777 (68.7)	<0.001
race.ethnicity (%)					<0.001
White	590 (52.2)	663 (58.7)	731 (64.6)	667 (59.0)
Black	172 (15.2)	91 (8.1)	74 (6.5)	103 (9.1)
Hispanic	242 (21.4)	223 (19.8)	202 (17.9)	222 (19.6)
Asian	13 (1.1)	36 (3.2)	27 (2.4)	27 (2.4)
Other	114 (10.1)	116 (10.3)	97 (8.6)	112 (9.9)
high.educ (%)					<0.001
<= 12 grades	50 (4.4)	42 (3.7)	39 (3.4)	45 (4.0)
HS Degree	390 (34.5)	273 (24.2)	242 (21.4)	261 (23.1)
College Degree	182 (16.1)	150 (13.3)	116 (10.3)	109 (9.6)
Bachelor	289 (25.5)	370 (32.7)	366 (32.4)	363 (32.1)
Higher	220 (19.4)	293 (25.9)	368 (32.5)	351 (31.0)
_miss	1 (0.1)	2 (0.2)	0 (0.0)	2 (0.2)
married = yes (%)	719 (63.5)	816 (72.2)	873 (77.2)	797 (70.5)	<0.001
HHInc (%)					<0.001
[<50K]	362 (32.0)	248 (21.9)	195 (17.2)	238 (21.0)
[50K - 100K]	332 (29.3)	333 (29.5)	307 (27.1)	291 (25.7)
[>100K]	345 (30.5)	457 (40.4)	547 (48.4)	512 (45.3)
_miss	93 (8.2)	92 (8.1)	82 (7.3)	90 (8.0)
anthro_1_height_in (mean (sd))	55.38 (3.28)	55.47 (3.18)	55.42 (3.04)	55.46 (3.13)	0.893
anthroweight1lb (mean (sd))	84.83 (24.05)	82.18 (22.88)	79.93 (21.37)	82.28 (22.75)	<0.001
demo_prnt_age_v2 (mean (sd))	39.73 (7.18)	40.35 (6.62)	41.02 (6.37)	40.72 (6.87)	<0.001
week_y_tvmovie (mean (sd))	14.07 (8.43)	9.22 (6.51)	6.56 (5.18)	5.58 (5.40)	<0.001
week_y_video (mean (sd))	12.46 (9.42)	5.51 (6.07)	3.26 (4.53)	3.69 (5.71)	<0.001
week_y_games (mean (sd))	14.68 (8.71)	6.14 (5.34)	3.68 (3.97)	3.32 (5.22)	<0.001
week_y_text (mean (sd))	0.94 (2.12)	0.94 (2.18)	1.02 (2.31)	2.98 (5.67)	<0.001
week_y_socnet (mean (sd))	0.31 (1.05)	0.35 (1.15)	0.36 (1.26)	1.80 (4.61)	<0.001
week_y_chat (mean (sd))	0.64 (1.72)	0.67 (1.81)	0.81 (1.77)	2.49 (4.92)	<0.001
kidsitotal (mean (sd))	0.39 (1.21)	0.27 (0.95)	0.20 (0.81)	0.22 (0.78)	<0.001
parsitotal (mean (sd))	0.01 (0.24)	0.03 (0.26)	0.03 (0.28)	0.03 (0.31)	0.611
KidsSIyes = 1 (%)	193 (17.3)	155 (13.9)	118 (10.6)	128 (11.4)	<0.001
ParSIyes = 1 (%)	8 (0.7)	15 (1.3)	13 (1.2)	13 (1.2)	0.534

Quartiles of Robust GFA5

	[-4.54,-0.612]	(-0.612,0.0165]	(0.0165,0.623]	(0.623,3.61]	p	test
n	1131	1132	1130	1131
age (mean (sd))	120.47 (7.23)	119.98 (7.46)	119.90 (7.40)	119.73 (7.25)	0.096
female = yes (%)	514 (45.4)	573 (50.6)	522 (46.2)	543 (48.0)	0.066
race.ethnicity (%)					<0.001
White	613 (54.2)	630 (55.8)	671 (59.4)	737 (65.2)
Black	131 (11.6)	121 (10.7)	104 (9.2)	84 (7.4)
Hispanic	249 (22.0)	251 (22.2)	206 (18.2)	183 (16.2)
Asian	27 (2.4)	27 (2.4)	32 (2.8)	17 (1.5)
Other	111 (9.8)	101 (8.9)	117 (10.4)	110 (9.7)
high.educ (%)					<0.001
<= 12 grades	43 (3.8)	67 (5.9)	36 (3.2)	30 (2.7)
HS Degree	330 (29.2)	299 (26.4)	302 (26.7)	235 (20.8)
College Degree	138 (12.2)	136 (12.0)	155 (13.7)	128 (11.3)
Bachelor	332 (29.4)	341 (30.1)	349 (30.9)	366 (32.4)
Higher	287 (25.4)	287 (25.4)	286 (25.3)	372 (32.9)
_miss	1 (0.1)	2 (0.2)	2 (0.2)	0 (0.0)
married = yes (%)	772 (68.3)	776 (68.6)	808 (71.5)	849 (75.1)	0.001
HHInc (%)					<0.001
[<50K]	295 (26.1)	288 (25.4)	244 (21.6)	216 (19.1)
[50K - 100K]	311 (27.5)	310 (27.4)	338 (29.9)	304 (26.9)
[>100K]	430 (38.0)	434 (38.3)	456 (40.4)	541 (47.8)
_miss	95 (8.4)	100 (8.8)	92 (8.1)	70 (6.2)
anthro_1_height_in (mean (sd))	55.54 (3.17)	55.52 (3.17)	55.31 (3.15)	55.35 (3.15)	0.187
anthroweight1lb (mean (sd))	83.37 (24.11)	83.63 (23.39)	81.74 (22.30)	80.47 (21.34)	0.003
demo_prnt_age_v2 (mean (sd))	39.87 (6.91)	40.12 (6.75)	40.76 (7.15)	41.06 (6.23)	<0.001
week_y_tvmovie (mean (sd))	8.97 (7.23)	9.00 (7.42)	8.91 (7.28)	8.54 (7.25)	0.417
week_y_video (mean (sd))	6.48 (7.74)	6.40 (7.82)	6.45 (7.61)	5.60 (7.35)	0.017
week_y_games (mean (sd))	7.40 (7.81)	6.80 (7.46)	7.22 (7.77)	6.41 (7.37)	0.009
week_y_text (mean (sd))	1.65 (3.92)	1.61 (3.84)	1.35 (3.20)	1.29 (3.06)	0.029
week_y_socnet (mean (sd))	0.70 (2.46)	0.94 (3.30)	0.60 (2.16)	0.58 (2.28)	0.003
week_y_chat (mean (sd))	1.15 (3.02)	1.24 (3.19)	1.11 (2.78)	1.12 (3.00)	0.729
kidsitotal (mean (sd))	0.27 (0.94)	0.28 (1.00)	0.25 (0.93)	0.28 (0.95)	0.780
parsitotal (mean (sd))	0.03 (0.22)	0.03 (0.40)	0.02 (0.23)	0.02 (0.19)	0.414
KidsSIyes = 1 (%)	145 (13.0)	155 (13.8)	144 (12.8)	150 (13.5)	0.893
ParSIyes = 1 (%)	17 (1.5)	11 (1.0)	11 (1.0)	10 (0.9)	0.447

Quartiles of Robust GFA6

	[-5.54,-0.462]	(-0.462,-0.051]	(-0.051,0.485]	(0.485,3.85]	p	test
n	1131	1132	1130	1131
age (mean (sd))	119.43 (7.51)	120.16 (7.27)	119.99 (7.21)	120.51 (7.32)	0.005
female = yes (%)	508 (44.9)	558 (49.3)	560 (49.6)	526 (46.5)	0.077
race.ethnicity (%)					0.002
White	658 (58.2)	666 (58.9)	652 (57.7)	675 (59.7)
Black	131 (11.6)	116 (10.3)	104 (9.2)	89 (7.9)
Hispanic	192 (17.0)	220 (19.5)	240 (21.2)	237 (21.0)
Asian	22 (1.9)	37 (3.3)	28 (2.5)	16 (1.4)
Other	127 (11.2)	92 (8.1)	106 (9.4)	114 (10.1)
high.educ (%)					0.109
<= 12 grades	48 (4.2)	40 (3.5)	36 (3.2)	52 (4.6)
HS Degree	305 (27.0)	286 (25.3)	303 (26.8)	272 (24.0)
College Degree	121 (10.7)	138 (12.2)	132 (11.7)	166 (14.7)
Bachelor	348 (30.8)	347 (30.7)	335 (29.6)	358 (31.7)
Higher	308 (27.2)	318 (28.1)	324 (28.7)	282 (24.9)
_miss	1 (0.1)	3 (0.3)	0 (0.0)	1 (0.1)
married = yes (%)	762 (67.4)	821 (72.5)	831 (73.5)	791 (69.9)	0.006
HHInc (%)					0.040
[<50K]	274 (24.2)	238 (21.0)	257 (22.7)	274 (24.2)
[50K - 100K]	292 (25.8)	306 (27.0)	330 (29.2)	335 (29.6)
[>100K]	458 (40.5)	499 (44.1)	468 (41.4)	436 (38.5)
_miss	107 (9.5)	89 (7.9)	75 (6.6)	86 (7.6)
anthro_1_height_in (mean (sd))	55.28 (3.19)	55.54 (3.00)	55.44 (3.22)	55.46 (3.22)	0.257
anthroweight1lb (mean (sd))	80.87 (21.13)	82.02 (21.97)	83.24 (24.42)	83.08 (23.65)	0.050
demo_prnt_age_v2 (mean (sd))	40.53 (6.96)	40.62 (7.02)	40.55 (6.66)	40.12 (6.48)	0.285
week_y_tvmovie (mean (sd))	9.19 (7.73)	8.33 (6.80)	8.97 (7.27)	8.94 (7.33)	0.034
week_y_video (mean (sd))	6.27 (7.74)	5.51 (7.04)	5.89 (7.30)	7.25 (8.31)	<0.001
week_y_games (mean (sd))	7.47 (8.10)	6.39 (7.08)	6.69 (7.37)	7.28 (7.82)	0.002
week_y_text (mean (sd))	1.59 (3.90)	1.48 (3.55)	1.43 (3.39)	1.39 (3.23)	0.588
week_y_socnet (mean (sd))	0.87 (3.19)	0.62 (2.12)	0.64 (2.39)	0.69 (2.54)	0.095
week_y_chat (mean (sd))	1.28 (3.63)	1.05 (2.31)	1.13 (2.96)	1.15 (2.95)	0.342
kidsitotal (mean (sd))	0.34 (1.08)	0.19 (0.75)	0.23 (0.86)	0.32 (1.08)	<0.001
parsitotal (mean (sd))	0.04 (0.38)	0.01 (0.12)	0.01 (0.17)	0.04 (0.34)	0.003
KidsSIyes = 1 (%)	182 (16.3)	113 (10.1)	134 (12.0)	165 (14.7)	<0.001
ParSIyes = 1 (%)	14 (1.3)	7 (0.6)	4 (0.4)	24 (2.2)	<0.001

Quartiles of Robust GFA7

	[-4.73,-0.523]	(-0.523,-0.101]	(-0.101,0.532]	(0.532,3.45]	p	test
n	1132	1130	1132	1130
age (mean (sd))	119.91 (7.35)	120.19 (7.36)	120.09 (7.28)	119.90 (7.36)	0.742
female = yes (%)	450 (39.8)	553 (48.9)	549 (48.5)	600 (53.1)	<0.001
race.ethnicity (%)					0.874
White	669 (59.2)	666 (58.9)	670 (59.2)	646 (57.2)
Black	99 (8.8)	109 (9.6)	110 (9.7)	122 (10.8)
Hispanic	229 (20.2)	221 (19.6)	213 (18.8)	226 (20.0)
Asian	30 (2.7)	28 (2.5)	25 (2.2)	20 (1.8)
Other	104 (9.2)	106 (9.4)	113 (10.0)	116 (10.3)
high.educ (%)					0.013
<= 12 grades	55 (4.9)	43 (3.8)	32 (2.8)	46 (4.1)
HS Degree	287 (25.4)	256 (22.7)	295 (26.1)	328 (29.0)
College Degree	121 (10.7)	146 (12.9)	142 (12.5)	148 (13.1)
Bachelor	348 (30.7)	373 (33.0)	332 (29.3)	335 (29.6)
Higher	319 (28.2)	312 (27.6)	329 (29.1)	272 (24.1)
_miss	2 (0.2)	0 (0.0)	2 (0.2)	1 (0.1)
married = yes (%)	833 (73.6)	837 (74.1)	789 (69.7)	746 (66.0)	<0.001
HHInc (%)					0.007
[<50K]	260 (23.0)	230 (20.4)	269 (23.8)	284 (25.1)
[50K - 100K]	282 (24.9)	326 (28.8)	327 (28.9)	328 (29.0)
[>100K]	479 (42.3)	486 (43.0)	463 (40.9)	433 (38.3)
_miss	111 (9.8)	88 (7.8)	73 (6.4)	85 (7.5)
anthro_1_height_in (mean (sd))	55.28 (3.27)	55.50 (3.10)	55.47 (3.11)	55.47 (3.16)	0.314
anthroweight1lb (mean (sd))	80.91 (22.70)	82.82 (23.80)	82.40 (23.09)	83.08 (21.70)	0.108
demo_prnt_age_v2 (mean (sd))	40.39 (6.88)	40.70 (6.84)	40.47 (6.45)	40.27 (6.97)	0.493
week_y_tvmovie (mean (sd))	8.28 (7.16)	8.57 (6.97)	9.02 (7.48)	9.56 (7.51)	<0.001
week_y_video (mean (sd))	5.65 (7.38)	5.80 (7.12)	6.22 (7.57)	7.25 (8.33)	<0.001
week_y_games (mean (sd))	7.16 (7.73)	6.74 (7.38)	6.98 (7.67)	6.94 (7.66)	0.627
week_y_text (mean (sd))	1.77 (4.49)	1.34 (3.03)	1.41 (3.32)	1.38 (3.06)	0.013
week_y_socnet (mean (sd))	0.80 (3.28)	0.62 (2.10)	0.68 (2.37)	0.72 (2.46)	0.378
week_y_chat (mean (sd))	1.41 (3.87)	1.01 (2.27)	1.07 (2.72)	1.11 (2.89)	0.007
kidsitotal (mean (sd))	0.35 (1.08)	0.19 (0.77)	0.21 (0.80)	0.33 (1.11)	<0.001
parsitotal (mean (sd))	0.04 (0.38)	0.01 (0.14)	0.01 (0.19)	0.03 (0.32)	0.020
KidsSIyes = 1 (%)	175 (15.6)	116 (10.4)	125 (11.2)	178 (15.9)	<0.001
ParSIyes = 1 (%)	20 (1.8)	6 (0.5)	7 (0.6)	16 (1.4)	0.010

Quartiles of Robust GFA8

	[-2.86,-0.458]	(-0.458,0.0765]	(0.0765,0.491]	(0.491,3.35]	p	test
n	1131	1133	1129	1131
age (mean (sd))	119.89 (7.37)	120.01 (7.43)	120.13 (7.41)	120.06 (7.14)	0.892
female = yes (%)	587 (51.9)	614 (54.2)	554 (49.1)	397 (35.1)	<0.001
race.ethnicity (%)					<0.001
White	659 (58.3)	700 (61.8)	707 (62.6)	585 (51.8)
Black	145 (12.8)	97 (8.6)	79 (7.0)	119 (10.5)
Hispanic	202 (17.9)	217 (19.2)	183 (16.2)	287 (25.4)
Asian	16 (1.4)	22 (1.9)	42 (3.7)	23 (2.0)
Other	109 (9.6)	97 (8.6)	118 (10.5)	115 (10.2)
high.educ (%)					<0.001
<= 12 grades	46 (4.1)	38 (3.4)	36 (3.2)	56 (5.0)
HS Degree	323 (28.6)	235 (20.7)	247 (21.9)	361 (31.9)
College Degree	149 (13.2)	139 (12.3)	118 (10.5)	151 (13.4)
Bachelor	345 (30.5)	366 (32.3)	367 (32.5)	310 (27.4)
Higher	266 (23.5)	354 (31.2)	360 (31.9)	252 (22.3)
_miss	2 (0.2)	1 (0.1)	1 (0.1)	1 (0.1)
married = yes (%)	752 (66.5)	843 (74.4)	852 (75.5)	758 (67.0)	<0.001
HHInc (%)					<0.001
[<50K]	296 (26.2)	227 (20.0)	221 (19.6)	299 (26.4)
[50K - 100K]	346 (30.6)	295 (26.0)	294 (26.0)	328 (29.0)
[>100K]	403 (35.6)	521 (46.0)	539 (47.7)	398 (35.2)
_miss	86 (7.6)	90 (7.9)	75 (6.6)	106 (9.4)
anthro_1_height_in (mean (sd))	55.51 (3.18)	55.44 (3.19)	55.49 (3.03)	55.29 (3.23)	0.343
anthroweight1lb (mean (sd))	84.87 (24.54)	80.54 (20.96)	80.34 (22.05)	83.46 (23.35)	<0.001
demo_prnt_age_v2 (mean (sd))	39.86 (7.07)	40.80 (6.34)	41.00 (6.79)	40.16 (6.87)	<0.001
week_y_tvmovie (mean (sd))	17.25 (6.86)	8.49 (5.05)	5.41 (4.60)	4.27 (3.95)	<0.001
week_y_video (mean (sd))	5.07 (6.93)	4.58 (6.42)	4.53 (5.86)	10.75 (9.08)	<0.001
week_y_games (mean (sd))	7.41 (7.90)	5.64 (6.55)	5.26 (6.29)	9.51 (8.69)	<0.001
week_y_text (mean (sd))	2.10 (4.90)	1.28 (3.03)	1.09 (2.62)	1.41 (3.02)	<0.001
week_y_socnet (mean (sd))	0.95 (3.39)	0.65 (2.52)	0.52 (2.04)	0.71 (2.20)	0.001
week_y_chat (mean (sd))	1.12 (3.20)	1.01 (2.27)	0.95 (2.42)	1.54 (3.82)	<0.001
kidsitotal (mean (sd))	0.29 (1.00)	0.26 (0.96)	0.27 (1.02)	0.26 (0.83)	0.899
parsitotal (mean (sd))	0.03 (0.26)	0.04 (0.37)	0.01 (0.16)	0.02 (0.27)	0.176
KidsSIyes = 1 (%)	153 (13.7)	145 (12.9)	137 (12.3)	159 (14.2)	0.549
ParSIyes = 1 (%)	15 (1.3)	16 (1.4)	6 (0.5)	12 (1.1)	0.168

# by z-scores of the factors < -.5, -.5 - .5, > .5
mysdGFA <- paste0("sdGFA",rep(1:8))
for(i in 1:length(mysdGFA)){
currdata[,mysdGFA[i]] <- cut(currdata[,paste0("SMA_RGFA",i)], breaks = rep(-2.5:3.5))
}

mynames <- names(currdata)

# General Sample Characteristics:
# Need to rename the demographic variables:
demovars <- c("age","female","race.ethnicity","high.educ","married",
              "HHInc","anthro_1_height_in","anthroweight1lb","demo_prnt_age_v2")

# Additional Screen variables - youth reported hours a week:
addscreenvars <- mynames[grep("week_y_",mynames)]

listvars <- c(demovars,addscreenvars,"kidsitotal","parsitotal","KidsSIyes","ParSIyes")

tabledata <- currdata[,c(mysdGFA,listvars)]

catVars = c("female","race.ethnicity","married","high.educ","HHInc","KidsSIyes","ParSIyes")
# length(myqGFA)
for(i in 1: length(mysdGFA)){
mytable1 <- CreateTableOne(vars = listvars,data=tabledata,factorVars=catVars,strata=c(mysdGFA[i]))

# Make the Table look nice:
tabAsStringMatrix <- print(mytable1, printToggle = FALSE, noSpaces = TRUE) 
print(kable(tabAsStringMatrix, "html",caption = paste0("Standard Deviations of Robust GFA",i)) %>% kable_styling(bootstrap_options = c( "hover", "condensed", "responsive"),font_size = 11))

}

Standard Deviations of Robust GFA1

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	1	1529	1996	665	228	81
age (mean (sd))	108.00 (NA)	120.40 (7.20)	119.82 (7.38)	119.82 (7.51)	119.78 (7.12)	120.57 (7.77)	NA
female = yes (%)	1 (100.0)	797 (52.1)	968 (48.5)	266 (40.0)	90 (39.5)	23 (28.4)	<0.001
race.ethnicity (%)							<0.001
White	0 (0.0)	855 (55.9)	1234 (61.9)	377 (56.8)	133 (58.3)	44 (54.3)
Black	1 (100.0)	180 (11.8)	157 (7.9)	58 (8.7)	28 (12.3)	11 (13.6)
Hispanic	0 (0.0)	330 (21.6)	356 (17.8)	143 (21.5)	40 (17.5)	13 (16.0)
Asian	0 (0.0)	37 (2.4)	53 (2.7)	11 (1.7)	2 (0.9)	0 (0.0)
Other	0 (0.0)	127 (8.3)	195 (9.8)	75 (11.3)	25 (11.0)	13 (16.0)
high.educ (%)							<0.001
<= 12 grades	0 (0.0)	62 (4.1)	56 (2.8)	39 (5.9)	10 (4.4)	7 (8.6)
HS Degree	0 (0.0)	387 (25.3)	479 (24.0)	197 (29.6)	62 (27.2)	29 (35.8)
College Degree	0 (0.0)	202 (13.2)	215 (10.8)	85 (12.8)	35 (15.4)	16 (19.8)
Bachelor	1 (100.0)	460 (30.1)	645 (32.3)	186 (28.0)	74 (32.5)	18 (22.2)
Higher	0 (0.0)	417 (27.3)	600 (30.1)	156 (23.5)	47 (20.6)	10 (12.3)
_miss	0 (0.0)	1 (0.1)	1 (0.1)	2 (0.3)	0 (0.0)	1 (1.2)
married = yes (%)	1 (100.0)	1126 (73.6)	1439 (72.1)	443 (66.6)	145 (63.6)	41 (50.6)	<0.001
HHInc (%)							<0.001
[<50K]	0 (0.0)	329 (21.5)	405 (20.3)	197 (29.6)	72 (31.6)	29 (35.8)
[50K - 100K]	1 (100.0)	418 (27.3)	568 (28.5)	185 (27.8)	65 (28.5)	20 (24.7)
[>100K]	0 (0.0)	646 (42.2)	888 (44.5)	233 (35.0)	72 (31.6)	19 (23.5)
_miss	0 (0.0)	136 (8.9)	135 (6.8)	50 (7.5)	19 (8.3)	13 (16.0)
anthro_1_height_in (mean (sd))	55.50 (NA)	55.52 (3.11)	55.41 (3.17)	55.35 (3.23)	55.27 (3.18)	55.23 (3.19)	NA
anthroweight1lb (mean (sd))	106.50 (NA)	82.82 (23.02)	81.24 (21.94)	84.16 (24.47)	82.37 (24.29)	83.08 (23.72)	NA
demo_prnt_age_v2 (mean (sd))	34.00 (NA)	40.58 (6.81)	40.63 (6.56)	40.11 (6.85)	40.01 (7.24)	38.12 (7.89)	NA
week_y_tvmovie (mean (sd))	10.50 (NA)	8.64 (7.06)	8.44 (7.04)	9.63 (7.66)	10.59 (8.77)	10.80 (7.83)	NA
week_y_video (mean (sd))	4.50 (NA)	5.50 (6.85)	6.17 (7.66)	7.02 (7.99)	7.70 (9.24)	9.48 (9.85)	NA
week_y_games (mean (sd))	2.50 (NA)	6.24 (7.08)	6.78 (7.40)	7.80 (8.09)	8.90 (9.11)	10.92 (9.65)	NA
week_y_text (mean (sd))	24.00 (NA)	1.55 (3.56)	1.26 (3.13)	1.49 (3.49)	2.43 (5.39)	1.77 (4.07)	NA
week_y_socnet (mean (sd))	0.00 (NA)	0.65 (2.22)	0.60 (2.27)	0.86 (3.09)	1.46 (4.81)	0.39 (1.38)	NA
week_y_chat (mean (sd))	11.00 (NA)	1.16 (2.87)	1.05 (2.74)	1.22 (3.19)	1.72 (4.74)	0.88 (3.03)	NA
kidsitotal (mean (sd))	0.00 (NA)	0.14 (0.58)	0.23 (0.88)	0.43 (1.21)	0.61 (1.52)	0.96 (2.04)	NA
parsitotal (mean (sd))	0.00 (NA)	0.00 (0.06)	0.01 (0.14)	0.04 (0.26)	0.17 (0.80)	0.23 (0.92)	NA
KidsSIyes = 1 (%)	0 (0.0)	128 (8.5)	232 (11.7)	129 (19.5)	63 (28.0)	30 (38.0)	<0.001
ParSIyes = 1 (%)	0 (0.0)	3 (0.2)	7 (0.4)	16 (2.4)	14 (6.3)	6 (7.6)	<0.001

Standard Deviations of Robust GFA2

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	236	1085	1806	1077	263	21
age (mean (sd))	120.83 (6.95)	120.91 (7.26)	119.86 (7.39)	119.47 (7.28)	118.49 (7.40)	118.43 (7.08)	<0.001
female = yes (%)	123 (52.1)	546 (50.3)	875 (48.4)	476 (44.2)	106 (40.3)	6 (28.6)	0.002
race.ethnicity (%)							<0.001
White	167 (70.8)	760 (70.0)	1121 (62.1)	506 (47.1)	72 (27.4)	2 (9.5)
Black	1 (0.4)	38 (3.5)	132 (7.3)	176 (16.4)	82 (31.2)	11 (52.4)
Hispanic	23 (9.7)	140 (12.9)	362 (20.0)	274 (25.5)	80 (30.4)	6 (28.6)
Asian	17 (7.2)	38 (3.5)	31 (1.7)	12 (1.1)	2 (0.8)	0 (0.0)
Other	28 (11.9)	109 (10.0)	160 (8.9)	107 (10.0)	27 (10.3)	2 (9.5)
high.educ (%)							<0.001
<= 12 grades	0 (0.0)	6 (0.6)	52 (2.9)	78 (7.2)	33 (12.5)	7 (33.3)
HS Degree	23 (9.7)	172 (15.9)	427 (23.6)	409 (38.0)	125 (47.5)	6 (28.6)
College Degree	12 (5.1)	101 (9.3)	232 (12.8)	172 (16.0)	34 (12.9)	5 (23.8)
Bachelor	69 (29.2)	388 (35.8)	620 (34.3)	259 (24.0)	38 (14.4)	0 (0.0)
Higher	132 (55.9)	417 (38.4)	474 (26.2)	157 (14.6)	33 (12.5)	2 (9.5)
_miss	0 (0.0)	1 (0.1)	1 (0.1)	2 (0.2)	0 (0.0)	1 (4.8)
married = yes (%)	208 (88.1)	861 (79.4)	1329 (73.6)	642 (59.6)	124 (47.1)	8 (38.1)	<0.001
HHInc (%)							<0.001
[<50K]	15 (6.4)	136 (12.5)	371 (20.5)	370 (34.4)	136 (51.7)	12 (57.1)
[50K - 100K]	53 (22.5)	298 (27.5)	545 (30.2)	300 (27.9)	56 (21.3)	2 (9.5)
[>100K]	151 (64.0)	589 (54.3)	761 (42.1)	297 (27.6)	37 (14.1)	2 (9.5)
_miss	17 (7.2)	62 (5.7)	129 (7.1)	110 (10.2)	34 (12.9)	5 (23.8)
anthro_1_height_in (mean (sd))	55.85 (3.13)	55.66 (2.96)	55.39 (3.10)	55.32 (3.15)	54.81 (4.07)	55.48 (3.77)	0.001
anthroweight1lb (mean (sd))	79.70 (19.63)	79.96 (19.81)	81.70 (22.20)	84.94 (25.08)	87.33 (28.43)	90.98 (34.03)	<0.001
demo_prnt_age_v2 (mean (sd))	42.34 (5.04)	41.72 (6.16)	40.57 (6.62)	39.11 (7.18)	37.96 (8.22)	37.57 (7.06)	<0.001
week_y_tvmovie (mean (sd))	5.58 (5.04)	7.30 (6.04)	8.82 (7.00)	10.59 (8.17)	11.51 (8.79)	13.71 (11.10)	<0.001
week_y_video (mean (sd))	3.29 (4.83)	4.83 (6.52)	6.19 (7.49)	7.66 (8.47)	9.03 (8.99)	11.35 (10.03)	<0.001
week_y_games (mean (sd))	4.35 (5.05)	5.86 (6.78)	6.90 (7.37)	7.99 (8.21)	9.71 (9.82)	13.57 (10.90)	<0.001
week_y_text (mean (sd))	0.72 (1.21)	1.07 (2.64)	1.33 (3.05)	1.92 (4.43)	2.66 (5.35)	6.51 (8.49)	<0.001
week_y_socnet (mean (sd))	0.23 (0.76)	0.43 (1.61)	0.64 (2.32)	0.98 (3.22)	1.26 (3.87)	5.57 (9.52)	<0.001
week_y_chat (mean (sd))	0.56 (1.68)	0.75 (1.83)	1.00 (2.35)	1.53 (3.74)	2.33 (5.17)	8.23 (11.07)	<0.001
kidsitotal (mean (sd))	0.30 (1.21)	0.25 (0.88)	0.23 (0.85)	0.32 (1.09)	0.38 (1.12)	0.14 (0.48)	0.066
parsitotal (mean (sd))	0.02 (0.20)	0.01 (0.13)	0.02 (0.21)	0.04 (0.35)	0.07 (0.60)	0.00 (0.00)	0.015
KidsSIyes = 1 (%)	30 (13.0)	142 (13.2)	214 (12.0)	156 (14.7)	44 (16.9)	2 (9.5)	0.167
ParSIyes = 1 (%)	3 (1.3)	7 (0.7)	17 (1.0)	16 (1.5)	6 (2.3)	0 (0.0)	0.177

Standard Deviations of Robust GFA3

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	14	1316	2350	591	137	50
age (mean (sd))	115.00 (5.33)	119.04 (7.29)	120.41 (7.36)	120.63 (7.21)	119.84 (7.10)	120.90 (7.20)	<0.001
female = yes (%)	5 (35.7)	715 (54.3)	1085 (46.2)	219 (37.1)	63 (46.0)	31 (62.0)	<0.001
race.ethnicity (%)							<0.001
White	6 (42.9)	817 (62.1)	1434 (61.0)	315 (53.3)	47 (34.3)	17 (34.0)
Black	4 (28.6)	92 (7.0)	187 (8.0)	87 (14.7)	29 (21.2)	13 (26.0)
Hispanic	3 (21.4)	270 (20.5)	416 (17.7)	128 (21.7)	43 (31.4)	13 (26.0)
Asian	0 (0.0)	28 (2.1)	64 (2.7)	10 (1.7)	1 (0.7)	0 (0.0)
Other	1 (7.1)	108 (8.2)	248 (10.6)	51 (8.6)	17 (12.4)	7 (14.0)
high.educ (%)							<0.001
<= 12 grades	3 (21.4)	72 (5.5)	63 (2.7)	26 (4.4)	5 (3.6)	4 (8.0)
HS Degree	5 (35.7)	294 (22.3)	558 (23.7)	195 (33.0)	59 (43.1)	23 (46.0)
College Degree	2 (14.3)	145 (11.0)	283 (12.0)	94 (15.9)	19 (13.9)	5 (10.0)
Bachelor	2 (14.3)	419 (31.8)	737 (31.4)	160 (27.1)	42 (30.7)	13 (26.0)
Higher	2 (14.3)	384 (29.2)	708 (30.1)	115 (19.5)	12 (8.8)	5 (10.0)
_miss	0 (0.0)	2 (0.2)	1 (0.0)	1 (0.2)	0 (0.0)	0 (0.0)
married = yes (%)	6 (42.9)	971 (73.8)	1729 (73.6)	367 (62.1)	82 (59.9)	19 (38.0)	<0.001
HHInc (%)							<0.001
[<50K]	8 (57.1)	280 (21.3)	465 (19.8)	185 (31.3)	51 (37.2)	25 (50.0)
[50K - 100K]	4 (28.6)	352 (26.7)	664 (28.3)	172 (29.1)	45 (32.8)	10 (20.0)
[>100K]	1 (7.1)	570 (43.3)	1051 (44.7)	185 (31.3)	31 (22.6)	9 (18.0)
_miss	1 (7.1)	114 (8.7)	170 (7.2)	49 (8.3)	10 (7.3)	6 (12.0)
anthro_1_height_in (mean (sd))	54.84 (4.94)	55.16 (3.11)	55.46 (3.07)	55.64 (3.42)	55.84 (3.03)	56.59 (2.67)	<0.001
anthroweight1lb (mean (sd))	81.16 (27.18)	80.30 (22.21)	81.60 (21.54)	85.34 (25.02)	89.40 (24.84)	92.14 (27.97)	<0.001
demo_prnt_age_v2 (mean (sd))	38.36 (7.69)	40.65 (6.79)	40.75 (6.60)	39.65 (7.21)	38.43 (7.07)	40.26 (7.11)	<0.001
week_y_tvmovie (mean (sd))	2.73 (1.92)	5.73 (5.00)	8.72 (6.72)	13.35 (8.43)	15.48 (9.24)	14.61 (8.80)	<0.001
week_y_video (mean (sd))	1.66 (2.20)	2.17 (2.92)	5.39 (6.10)	13.67 (9.23)	18.23 (9.17)	15.76 (10.13)	<0.001
week_y_games (mean (sd))	1.57 (1.98)	2.79 (3.11)	6.23 (6.14)	14.62 (9.23)	18.01 (9.05)	15.85 (9.65)	<0.001
week_y_text (mean (sd))	0.09 (0.33)	0.34 (0.86)	0.95 (1.47)	2.20 (2.98)	6.45 (6.77)	11.07 (8.61)	<0.001
week_y_socnet (mean (sd))	0.00 (0.00)	0.09 (0.41)	0.36 (0.96)	1.14 (2.11)	2.28 (3.46)	6.94 (7.80)	<0.001
week_y_chat (mean (sd))	0.04 (0.13)	0.27 (0.70)	0.76 (1.37)	1.80 (2.53)	4.59 (5.63)	7.65 (8.17)	<0.001
kidsitotal (mean (sd))	0.50 (0.85)	0.25 (0.91)	0.24 (0.86)	0.33 (1.08)	0.50 (1.71)	0.34 (1.04)	0.018
parsitotal (mean (sd))	0.14 (0.53)	0.03 (0.32)	0.01 (0.21)	0.04 (0.39)	0.02 (0.19)	0.00 (0.00)	0.109
KidsSIyes = 1 (%)	5 (35.7)	153 (11.8)	300 (12.9)	91 (15.6)	24 (17.8)	8 (16.0)	0.013
ParSIyes = 1 (%)	1 (7.1)	16 (1.2)	18 (0.8)	10 (1.7)	2 (1.5)	0 (0.0)	0.073

Standard Deviations of Robust GFA4

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	240	792	2239	1058	79	30
age (mean (sd))	119.40 (7.07)	119.98 (7.41)	119.96 (7.28)	120.21 (7.42)	121.00 (7.63)	120.63 (6.74)	0.519
female = yes (%)	49 (20.4)	238 (30.1)	1040 (46.4)	728 (68.8)	50 (63.3)	20 (66.7)	<0.001
race.ethnicity (%)							<0.001
White	123 (51.2)	418 (52.8)	1380 (61.7)	659 (62.3)	34 (43.0)	5 (16.7)
Black	46 (19.2)	106 (13.4)	167 (7.5)	69 (6.5)	16 (20.3)	8 (26.7)
Hispanic	41 (17.1)	179 (22.6)	423 (18.9)	194 (18.3)	19 (24.1)	12 (40.0)
Asian	5 (2.1)	7 (0.9)	60 (2.7)	30 (2.8)	1 (1.3)	0 (0.0)
Other	25 (10.4)	81 (10.2)	208 (9.3)	106 (10.0)	9 (11.4)	5 (16.7)
high.educ (%)							<0.001
<= 12 grades	14 (5.8)	30 (3.8)	82 (3.7)	36 (3.4)	4 (5.1)	4 (13.3)
HS Degree	90 (37.5)	268 (33.8)	512 (22.9)	221 (20.9)	25 (31.6)	15 (50.0)
College Degree	50 (20.8)	117 (14.8)	264 (11.8)	104 (9.8)	7 (8.9)	2 (6.7)
Bachelor	50 (20.8)	209 (26.4)	736 (32.9)	341 (32.2)	26 (32.9)	7 (23.3)
Higher	36 (15.0)	167 (21.1)	643 (28.7)	355 (33.6)	17 (21.5)	2 (6.7)
_miss	0 (0.0)	1 (0.1)	2 (0.1)	1 (0.1)	0 (0.0)	0 (0.0)
married = yes (%)	140 (58.3)	519 (65.5)	1668 (74.5)	771 (72.9)	45 (57.0)	15 (50.0)	<0.001
HHInc (%)							<0.001
[<50K]	83 (34.6)	238 (30.1)	446 (19.9)	193 (18.2)	26 (32.9)	17 (56.7)
[50K - 100K]	73 (30.4)	231 (29.2)	637 (28.5)	277 (26.2)	15 (19.0)	7 (23.3)
[>100K]	59 (24.6)	262 (33.1)	985 (44.0)	503 (47.5)	31 (39.2)	5 (16.7)
_miss	25 (10.4)	61 (7.7)	171 (7.6)	85 (8.0)	7 (8.9)	1 (3.3)
anthro_1_height_in (mean (sd))	54.99 (3.56)	55.50 (3.17)	55.48 (3.10)	55.34 (3.14)	55.62 (2.93)	57.31 (3.60)	0.004
anthroweight1lb (mean (sd))	84.11 (23.81)	85.11 (23.78)	81.23 (22.27)	81.16 (21.93)	84.81 (23.10)	97.73 (31.00)	<0.001
demo_prnt_age_v2 (mean (sd))	38.75 (6.94)	40.12 (7.32)	40.66 (6.53)	40.99 (6.71)	39.14 (7.56)	38.20 (7.03)	<0.001
week_y_tvmovie (mean (sd))	16.80 (8.28)	12.80 (8.05)	8.03 (6.08)	4.81 (4.30)	8.75 (6.85)	11.97 (9.09)	<0.001
week_y_video (mean (sd))	16.37 (9.47)	10.61 (8.69)	4.50 (5.55)	2.71 (3.95)	10.10 (9.75)	9.66 (9.54)	<0.001
week_y_games (mean (sd))	20.97 (7.64)	12.25 (7.57)	4.99 (4.91)	2.58 (3.68)	7.74 (8.56)	9.76 (8.74)	<0.001
week_y_text (mean (sd))	0.73 (1.59)	1.03 (2.32)	0.97 (2.16)	1.67 (3.26)	7.96 (6.81)	14.16 (9.58)	<0.001
week_y_socnet (mean (sd))	0.29 (0.92)	0.31 (1.10)	0.35 (1.21)	0.73 (1.65)	5.36 (6.64)	9.57 (8.04)	<0.001
week_y_chat (mean (sd))	0.71 (1.54)	0.65 (1.85)	0.73 (1.75)	1.42 (2.74)	6.94 (7.03)	11.72 (9.89)	<0.001
kidsitotal (mean (sd))	0.48 (1.31)	0.36 (1.14)	0.24 (0.88)	0.18 (0.67)	0.57 (1.47)	0.31 (0.85)	<0.001
parsitotal (mean (sd))	0.01 (0.09)	0.02 (0.27)	0.03 (0.28)	0.03 (0.30)	0.00 (0.00)	0.03 (0.18)	0.790
KidsSIyes = 1 (%)	48 (20.3)	129 (16.5)	274 (12.4)	108 (10.3)	16 (20.3)	4 (13.8)	<0.001
ParSIyes = 1 (%)	2 (0.9)	5 (0.6)	28 (1.3)	11 (1.1)	0 (0.0)	1 (3.3)	0.474

Standard Deviations of Robust GFA5

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	263	977	1914	1020	268	31
age (mean (sd))	121.52 (6.92)	120.02 (7.29)	119.97 (7.52)	119.63 (7.25)	120.13 (7.02)	120.87 (6.08)	0.014
female = yes (%)	110 (41.8)	465 (47.6)	933 (48.7)	482 (47.3)	125 (46.6)	15 (48.4)	0.456
race.ethnicity (%)							<0.001
White	147 (55.9)	531 (54.4)	1098 (57.4)	615 (60.3)	204 (76.1)	26 (83.9)
Black	36 (13.7)	102 (10.4)	196 (10.3)	84 (8.2)	15 (5.6)	1 (3.2)
Hispanic	49 (18.6)	223 (22.8)	387 (20.2)	195 (19.1)	24 (9.0)	1 (3.2)
Asian	10 (3.8)	22 (2.3)	47 (2.5)	22 (2.2)	2 (0.7)	0 (0.0)
Other	21 (8.0)	99 (10.1)	184 (9.6)	104 (10.2)	23 (8.6)	3 (9.7)
high.educ (%)							<0.001
<= 12 grades	12 (4.6)	42 (4.3)	87 (4.5)	35 (3.4)	0 (0.0)	0 (0.0)
HS Degree	64 (24.3)	284 (29.1)	515 (26.9)	229 (22.5)	55 (20.5)	2 (6.5)
College Degree	30 (11.4)	127 (13.0)	235 (12.3)	138 (13.5)	17 (6.3)	3 (9.7)
Bachelor	82 (31.2)	289 (29.6)	577 (30.1)	322 (31.6)	91 (34.0)	12 (38.7)
Higher	74 (28.1)	235 (24.1)	496 (25.9)	296 (29.0)	105 (39.2)	14 (45.2)
_miss	1 (0.4)	0 (0.0)	4 (0.2)	0 (0.0)	0 (0.0)	0 (0.0)
married = yes (%)	177 (67.3)	669 (68.5)	1332 (69.6)	751 (73.6)	217 (81.0)	26 (83.9)	<0.001
HHInc (%)							<0.001
[<50K]	71 (27.0)	248 (25.4)	460 (24.0)	217 (21.3)	35 (13.1)	4 (12.9)
[50K - 100K]	61 (23.2)	280 (28.7)	544 (28.4)	284 (27.8)	68 (25.4)	5 (16.1)
[>100K]	112 (42.6)	360 (36.8)	755 (39.4)	445 (43.6)	153 (57.1)	21 (67.7)
_miss	19 (7.2)	89 (9.1)	155 (8.1)	74 (7.3)	12 (4.5)	1 (3.2)
anthro_1_height_in (mean (sd))	55.76 (3.15)	55.45 (3.21)	55.44 (3.15)	55.34 (3.11)	55.27 (3.39)	55.13 (2.57)	0.462
anthroweight1lb (mean (sd))	83.81 (24.56)	82.88 (23.82)	82.97 (22.88)	81.37 (22.29)	78.48 (19.76)	76.34 (19.13)	0.012
demo_prnt_age_v2 (mean (sd))	40.39 (6.77)	39.70 (6.94)	40.40 (6.92)	40.79 (6.50)	42.00 (6.13)	43.13 (5.51)	<0.001
week_y_tvmovie (mean (sd))	8.82 (7.12)	9.15 (7.33)	8.88 (7.33)	8.75 (7.31)	8.38 (7.28)	8.03 (6.90)	0.654
week_y_video (mean (sd))	5.91 (7.22)	6.54 (7.89)	6.43 (7.71)	5.86 (7.49)	5.28 (7.08)	5.90 (6.69)	0.075
week_y_games (mean (sd))	6.47 (7.06)	7.60 (7.93)	6.95 (7.69)	6.68 (7.46)	6.07 (6.91)	7.13 (6.83)	0.022
week_y_text (mean (sd))	1.29 (3.22)	1.77 (4.10)	1.49 (3.58)	1.28 (3.07)	1.09 (2.07)	1.75 (5.17)	0.016
week_y_socnet (mean (sd))	0.49 (1.57)	0.81 (2.81)	0.78 (2.84)	0.60 (2.27)	0.44 (1.93)	0.38 (1.15)	0.074
week_y_chat (mean (sd))	0.76 (1.50)	1.31 (3.37)	1.16 (2.99)	1.05 (2.73)	1.27 (3.43)	1.25 (4.15)	0.114
kidsitotal (mean (sd))	0.35 (1.17)	0.26 (0.95)	0.26 (0.91)	0.28 (1.02)	0.26 (0.78)	0.35 (0.98)	0.716
parsitotal (mean (sd))	0.03 (0.25)	0.02 (0.21)	0.03 (0.35)	0.01 (0.17)	0.00 (0.06)	0.10 (0.54)	0.323
KidsSIyes = 1 (%)	38 (14.7)	121 (12.6)	250 (13.2)	132 (13.1)	38 (14.5)	6 (19.4)	0.813
ParSIyes = 1 (%)	5 (2.0)	13 (1.4)	18 (1.0)	10 (1.0)	1 (0.4)	1 (3.2)	0.370

Standard Deviations of Robust GFA6

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	155	820	2382	884	189	30
age (mean (sd))	119.19 (7.80)	119.60 (7.43)	120.07 (7.26)	120.56 (7.32)	120.01 (7.44)	119.47 (7.50)	0.082
female = yes (%)	66 (42.6)	382 (46.6)	1174 (49.3)	421 (47.6)	81 (42.9)	11 (36.7)	0.178
race.ethnicity (%)							0.042
White	93 (60.0)	480 (58.6)	1386 (58.2)	534 (60.4)	106 (56.1)	23 (76.7)
Black	19 (12.3)	81 (9.9)	236 (9.9)	75 (8.5)	9 (4.8)	1 (3.3)
Hispanic	23 (14.8)	150 (18.3)	481 (20.2)	180 (20.4)	45 (23.8)	5 (16.7)
Asian	1 (0.6)	19 (2.3)	67 (2.8)	12 (1.4)	4 (2.1)	0 (0.0)
Other	19 (12.3)	89 (10.9)	211 (8.9)	83 (9.4)	25 (13.2)	1 (3.3)
high.educ (%)							0.047
<= 12 grades	13 (8.4)	28 (3.4)	79 (3.3)	43 (4.9)	9 (4.8)	0 (0.0)
HS Degree	49 (31.6)	215 (26.2)	618 (25.9)	206 (23.3)	49 (25.9)	7 (23.3)
College Degree	18 (11.6)	83 (10.1)	283 (11.9)	132 (14.9)	25 (13.2)	6 (20.0)
Bachelor	42 (27.1)	263 (32.1)	720 (30.2)	276 (31.2)	63 (33.3)	11 (36.7)
Higher	33 (21.3)	230 (28.0)	679 (28.5)	226 (25.6)	43 (22.8)	6 (20.0)
_miss	0 (0.0)	1 (0.1)	3 (0.1)	1 (0.1)	0 (0.0)	0 (0.0)
married = yes (%)	105 (67.7)	554 (67.6)	1739 (73.0)	628 (71.0)	128 (67.7)	15 (50.0)	0.004
HHInc (%)							0.003
[<50K]	49 (31.6)	185 (22.6)	521 (21.9)	195 (22.1)	58 (30.7)	10 (33.3)
[50K - 100K]	35 (22.6)	221 (27.0)	665 (27.9)	266 (30.1)	55 (29.1)	8 (26.7)
[>100K]	54 (34.8)	333 (40.6)	1026 (43.1)	358 (40.5)	58 (30.7)	10 (33.3)
_miss	17 (11.0)	81 (9.9)	170 (7.1)	65 (7.4)	18 (9.5)	2 (6.7)
anthro_1_height_in (mean (sd))	54.86 (3.27)	55.47 (3.14)	55.48 (3.12)	55.37 (3.21)	55.74 (3.21)	55.66 (3.41)	0.168
anthroweight1lb (mean (sd))	78.74 (19.32)	81.49 (21.68)	82.62 (23.11)	82.15 (22.90)	85.40 (23.50)	86.61 (28.15)	0.078
demo_prnt_age_v2 (mean (sd))	39.87 (7.24)	40.75 (6.94)	40.59 (6.82)	40.12 (6.47)	40.19 (6.66)	39.00 (6.11)	0.186
week_y_tvmovie (mean (sd))	10.73 (8.57)	8.77 (7.53)	8.65 (7.03)	8.81 (7.28)	9.22 (7.52)	10.38 (7.53)	0.016
week_y_video (mean (sd))	7.10 (8.82)	6.09 (7.63)	5.74 (7.17)	7.15 (8.32)	7.07 (7.92)	9.49 (9.35)	<0.001
week_y_games (mean (sd))	9.51 (9.23)	6.98 (7.83)	6.59 (7.26)	7.18 (7.85)	6.98 (7.24)	10.10 (8.72)	<0.001
week_y_text (mean (sd))	1.78 (4.54)	1.54 (3.77)	1.45 (3.45)	1.32 (3.33)	1.53 (2.75)	1.10 (1.88)	0.627
week_y_socnet (mean (sd))	0.94 (3.14)	0.78 (2.88)	0.65 (2.38)	0.63 (2.45)	0.73 (2.10)	1.01 (3.19)	0.540
week_y_chat (mean (sd))	1.36 (3.12)	1.26 (3.66)	1.08 (2.62)	1.11 (2.85)	1.30 (3.55)	1.00 (1.85)	0.572
kidsitotal (mean (sd))	0.59 (1.58)	0.27 (0.89)	0.21 (0.80)	0.26 (0.95)	0.51 (1.56)	0.57 (1.04)	<0.001
parsitotal (mean (sd))	0.03 (0.20)	0.03 (0.35)	0.01 (0.14)	0.03 (0.28)	0.03 (0.23)	0.30 (1.21)	<0.001
KidsSIyes = 1 (%)	37 (23.9)	114 (14.1)	259 (11.0)	113 (12.9)	39 (20.7)	9 (30.0)	<0.001
ParSIyes = 1 (%)	3 (2.0)	7 (0.9)	11 (0.5)	16 (1.8)	4 (2.1)	2 (6.7)	<0.001

Standard Deviations of Robust GFA7

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	130	1023	2160	918	224	24
age (mean (sd))	119.82 (7.52)	120.07 (7.34)	120.10 (7.31)	119.93 (7.34)	120.10 (7.55)	118.62 (7.49)	0.926
female = yes (%)	31 (23.8)	450 (44.0)	1046 (48.4)	484 (52.7)	123 (54.9)	11 (45.8)	<0.001
race.ethnicity (%)							0.204
White	85 (65.4)	596 (58.3)	1276 (59.1)	536 (58.4)	130 (58.0)	6 (25.0)
Black	12 (9.2)	93 (9.1)	206 (9.5)	90 (9.8)	30 (13.4)	5 (20.8)
Hispanic	20 (15.4)	211 (20.6)	417 (19.3)	184 (20.0)	39 (17.4)	8 (33.3)
Asian	3 (2.3)	28 (2.7)	51 (2.4)	15 (1.6)	6 (2.7)	0 (0.0)
Other	10 (7.7)	94 (9.2)	209 (9.7)	93 (10.1)	19 (8.5)	5 (20.8)
high.educ (%)							0.030
<= 12 grades	5 (3.8)	51 (5.0)	71 (3.3)	37 (4.0)	8 (3.6)	3 (12.5)
HS Degree	31 (23.8)	266 (26.0)	521 (24.1)	270 (29.4)	58 (25.9)	8 (33.3)
College Degree	18 (13.8)	102 (10.0)	276 (12.8)	120 (13.1)	27 (12.1)	5 (20.8)
Bachelor	41 (31.5)	309 (30.2)	674 (31.2)	269 (29.3)	75 (33.5)	4 (16.7)
Higher	34 (26.2)	294 (28.7)	616 (28.5)	221 (24.1)	56 (25.0)	4 (16.7)
_miss	1 (0.8)	1 (0.1)	2 (0.1)	1 (0.1)	0 (0.0)	0 (0.0)
married = yes (%)	98 (75.4)	753 (73.6)	1554 (71.9)	618 (67.3)	142 (63.4)	11 (45.8)	<0.001
HHInc (%)							0.007
[<50K]	33 (25.4)	227 (22.2)	474 (21.9)	224 (24.4)	56 (25.0)	12 (50.0)
[50K - 100K]	34 (26.2)	253 (24.7)	625 (28.9)	269 (29.3)	66 (29.5)	5 (20.8)
[>100K]	50 (38.5)	443 (43.3)	909 (42.1)	355 (38.7)	86 (38.4)	4 (16.7)
_miss	13 (10.0)	100 (9.8)	152 (7.0)	70 (7.6)	16 (7.1)	3 (12.5)
anthro_1_height_in (mean (sd))	55.26 (3.14)	55.29 (3.25)	55.49 (3.13)	55.42 (3.10)	55.76 (3.22)	54.60 (3.45)	0.199
anthroweight1lb (mean (sd))	84.14 (30.42)	80.58 (21.38)	82.71 (23.61)	82.62 (21.75)	84.37 (21.40)	84.55 (18.96)	0.090
demo_prnt_age_v2 (mean (sd))	40.02 (7.14)	40.51 (6.82)	40.57 (6.63)	40.04 (6.88)	41.36 (6.96)	38.74 (8.34)	0.075
week_y_tvmovie (mean (sd))	9.89 (8.60)	8.02 (6.81)	8.83 (7.27)	9.50 (7.55)	9.56 (7.34)	11.91 (7.93)	<0.001
week_y_video (mean (sd))	6.16 (7.90)	5.46 (7.16)	6.08 (7.42)	7.10 (8.24)	7.55 (8.46)	7.39 (8.60)	<0.001
week_y_games (mean (sd))	8.90 (8.40)	6.79 (7.44)	6.92 (7.58)	6.89 (7.63)	6.90 (7.67)	8.25 (9.27)	0.081
week_y_text (mean (sd))	1.39 (3.98)	1.78 (4.48)	1.39 (3.18)	1.27 (2.96)	1.71 (3.40)	2.62 (5.03)	0.009
week_y_socnet (mean (sd))	0.60 (2.67)	0.82 (3.33)	0.67 (2.28)	0.69 (2.36)	0.75 (2.72)	0.75 (1.89)	0.736
week_y_chat (mean (sd))	1.22 (3.40)	1.44 (3.90)	1.05 (2.55)	1.08 (2.91)	1.19 (2.72)	1.32 (2.02)	0.030
kidsitotal (mean (sd))	0.53 (1.39)	0.30 (1.00)	0.20 (0.79)	0.32 (1.06)	0.34 (1.24)	0.29 (0.62)	<0.001
parsitotal (mean (sd))	0.09 (0.44)	0.03 (0.30)	0.01 (0.17)	0.03 (0.33)	0.02 (0.24)	0.04 (0.20)	0.026
KidsSIyes = 1 (%)	29 (22.7)	142 (14.0)	229 (10.8)	142 (15.6)	34 (15.2)	5 (20.8)	<0.001
ParSIyes = 1 (%)	6 (4.7)	11 (1.1)	12 (0.6)	14 (1.5)	2 (0.9)	1 (4.2)	<0.001

Standard Deviations of Robust GFA8

	(-2.5,-1.5]	(-1.5,-0.5]	(-0.5,0.5]	(0.5,1.5]	(1.5,2.5]	(2.5,3.5]	p	test
n	237	801	2355	958	136	23
age (mean (sd))	119.68 (7.27)	119.98 (7.41)	120.05 (7.40)	120.19 (7.18)	119.12 (7.05)	119.78 (5.74)	0.673
female = yes (%)	145 (61.2)	383 (47.8)	1221 (51.8)	346 (36.1)	42 (30.9)	4 (17.4)	<0.001
race.ethnicity (%)							<0.001
White	116 (48.9)	485 (60.5)	1465 (62.2)	495 (51.7)	72 (53.3)	10 (43.5)
Black	39 (16.5)	94 (11.7)	189 (8.0)	94 (9.8)	20 (14.8)	4 (17.4)
Hispanic	44 (18.6)	144 (18.0)	413 (17.5)	249 (26.0)	27 (20.0)	8 (34.8)
Asian	2 (0.8)	14 (1.7)	64 (2.7)	22 (2.3)	1 (0.7)	0 (0.0)
Other	36 (15.2)	64 (8.0)	224 (9.5)	97 (10.1)	15 (11.1)	1 (4.3)
high.educ (%)							<0.001
<= 12 grades	17 (7.2)	27 (3.4)	76 (3.2)	47 (4.9)	7 (5.1)	1 (4.3)
HS Degree	89 (37.6)	208 (26.0)	507 (21.5)	294 (30.7)	55 (40.4)	8 (34.8)
College Degree	35 (14.8)	110 (13.7)	263 (11.2)	124 (12.9)	15 (11.0)	9 (39.1)
Bachelor	46 (19.4)	264 (33.0)	766 (32.5)	265 (27.7)	38 (27.9)	5 (21.7)
Higher	50 (21.1)	190 (23.7)	741 (31.5)	227 (23.7)	21 (15.4)	0 (0.0)
_miss	0 (0.0)	2 (0.2)	2 (0.1)	1 (0.1)	0 (0.0)	0 (0.0)
married = yes (%)	143 (60.3)	547 (68.3)	1756 (74.6)	653 (68.2)	84 (61.8)	12 (52.2)	<0.001
HHInc (%)							<0.001
[<50K]	75 (31.6)	198 (24.7)	469 (19.9)	240 (25.1)	42 (30.9)	12 (52.2)
[50K - 100K]	75 (31.6)	249 (31.1)	613 (26.0)	277 (28.9)	42 (30.9)	6 (26.1)
[>100K]	67 (28.3)	291 (36.3)	1105 (46.9)	351 (36.6)	39 (28.7)	2 (8.7)
_miss	20 (8.4)	63 (7.9)	168 (7.1)	90 (9.4)	13 (9.6)	3 (13.0)
anthro_1_height_in (mean (sd))	55.74 (3.22)	55.42 (3.10)	55.46 (3.14)	55.33 (3.25)	54.92 (2.69)	56.13 (4.54)	0.145
anthroweight1lb (mean (sd))	88.46 (26.17)	83.75 (24.24)	80.57 (21.53)	83.52 (23.58)	81.69 (20.79)	92.12 (28.59)	<0.001
demo_prnt_age_v2 (mean (sd))	39.34 (7.11)	39.78 (6.94)	40.92 (6.57)	40.27 (6.77)	39.64 (7.78)	37.83 (6.51)	<0.001
week_y_tvmovie (mean (sd))	23.96 (3.87)	15.62 (6.11)	7.11 (5.20)	4.50 (4.02)	3.03 (3.08)	1.58 (1.87)	<0.001
week_y_video (mean (sd))	3.03 (4.35)	5.62 (7.31)	4.62 (6.24)	9.14 (8.17)	19.70 (7.81)	27.70 (1.11)	<0.001
week_y_games (mean (sd))	8.01 (8.12)	7.37 (7.89)	5.46 (6.44)	8.60 (8.14)	14.71 (9.44)	18.76 (10.84)	<0.001
week_y_text (mean (sd))	2.63 (5.67)	1.94 (4.56)	1.19 (2.84)	1.37 (2.94)	1.41 (2.65)	3.58 (6.50)	<0.001
week_y_socnet (mean (sd))	1.06 (3.50)	0.91 (3.32)	0.58 (2.27)	0.68 (2.24)	0.82 (1.84)	1.54 (2.84)	0.002
week_y_chat (mean (sd))	1.03 (2.78)	1.13 (3.25)	0.98 (2.39)	1.34 (3.35)	2.39 (5.24)	5.26 (8.24)	<0.001
kidsitotal (mean (sd))	0.33 (0.95)	0.29 (1.06)	0.26 (0.97)	0.26 (0.84)	0.32 (0.79)	0.17 (0.65)	0.821
parsitotal (mean (sd))	0.02 (0.20)	0.03 (0.28)	0.02 (0.28)	0.02 (0.29)	0.00 (0.00)	0.04 (0.21)	0.918
KidsSIyes = 1 (%)	41 (17.5)	102 (12.9)	292 (12.5)	130 (13.7)	27 (20.0)	2 (8.7)	0.059
ParSIyes = 1 (%)	3 (1.3)	11 (1.4)	23 (1.0)	11 (1.2)	0 (0.0)	1 (4.3)	0.453

# Set up the variables:
covars <- c("age","female","race.ethnicity","high.educ","married","HHInc","demo_prnt_age_v2")
colabels <- c("Age","Female","Race: Black","Race: Hispanic","Race: Asian","Race: Other","Parental Education: HS","Parental Education: College","Parental Education: Bachelor","Parental Education: > Bachelor","Married","Household Income: 50-100K","Household Income: > 100K","Household Income: miss","Parental Age")

## Select the nesting variables: site and twin status
nestvars <- c("site_name","FamilyID")
# Independent variables:
indepvars <- paste0("SMA_RGFA",rep(1:8))
# dependent variables:
depvars <- c("kidsitotal","parsitotal")

GFAselectLabel <- paste0("SMA_RGFA",rep(1:8))

# Creating a data frame:
glmmadmbdata <- currdata[complete.cases(currdata[,c(covars,nestvars,depvars,indepvars)]),c(covars,nestvars,depvars,indepvars)]

# Visualizing missing data
aggr(currdata[,c(covars,nestvars,depvars,indepvars)],col = c("blue","orange"),sortVars=TRUE,prop = FALSE, numbers = TRUE, combined = TRUE, cex.lab = 0.5, cex.axis =0.5, cex.numbers =0.5)

## 
##  Variables sorted by number of missings: 
##          Variable Count
##        parsitotal    70
##        kidsitotal    53
##  demo_prnt_age_v2    43
##    race.ethnicity     2
##               age     0
##            female     0
##         high.educ     0
##           married     0
##             HHInc     0
##         site_name     0
##          FamilyID     0
##         SMA_RGFA1     0
##         SMA_RGFA2     0
##         SMA_RGFA3     0
##         SMA_RGFA4     0
##         SMA_RGFA5     0
##         SMA_RGFA6     0
##         SMA_RGFA7     0
##         SMA_RGFA8     0

mydepvar <- depvars[1]

# Poisson Model with standard zero inflation:

kidsi_cov_ZIP <- myglmmTMBall(mydepvar,"null",covars,nestvars,glmmadmbdata,"poisson",1)
summary(kidsi_cov_ZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ age + female + race.ethnicity + high.educ + married +  
##     HHInc + demo_prnt_age_v2 + (1 | site_name)
## Zero inflation:              ~1
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5295.7   5417.1  -2628.9   5257.7     4376 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.03835  0.1958  
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.287713   0.666468  -0.432  0.66596    
## age                      0.008809   0.005010   1.758  0.07870 .  
## femaleyes               -0.153409   0.075186  -2.040  0.04131 *  
## race.ethnicityBlack     -0.092226   0.140903  -0.655  0.51277    
## race.ethnicityHispanic  -0.601359   0.123956  -4.851 1.23e-06 ***
## race.ethnicityAsian      0.557539   0.176700   3.155  0.00160 ** 
## race.ethnicityOther     -0.106349   0.114009  -0.933  0.35092    
## high.educHS Degree       0.322970   0.268911   1.201  0.22974    
## high.educCollege Degree  0.524344   0.279870   1.874  0.06100 .  
## high.educBachelor        0.427334   0.275599   1.551  0.12101    
## high.educHigher          0.267492   0.281748   0.949  0.34242    
## high.educ_miss           0.535192   0.836874   0.640  0.52249    
## marriedyes              -0.228173   0.087056  -2.621  0.00877 ** 
## HHInc[50K - 100K]       -0.243370   0.104648  -2.326  0.02004 *  
## HHInc[>100K]            -0.237761   0.113765  -2.090  0.03662 *  
## HHInc_miss              -0.433070   0.169179  -2.560  0.01047 *  
## demo_prnt_age_v2        -0.005046   0.005924  -0.852  0.39437    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.54697    0.05274   29.33   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(kidsi_cov_ZIP,title = "ZIP Covariates: Youth Suicide Items, standard zero inflation")

kidsi_sma_ZIP <- myglmmTMBall(mydepvar,indepvars,covars,nestvars,glmmadmbdata,"poisson",1)
summary(kidsi_sma_ZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 +  
##     SMA_RGFA5 + SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female +  
##     race.ethnicity + high.educ + married + HHInc + demo_prnt_age_v2 +  
##     (1 | site_name)
## Zero inflation:              ~1
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5183.7   5356.2  -2564.8   5129.7     4368 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.0649   0.2548  
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.631630   0.664875  -0.950   0.3421    
## SMA_RGFA1                0.270817   0.030413   8.905  < 2e-16 ***
## SMA_RGFA2                0.179537   0.039740   4.518 6.25e-06 ***
## SMA_RGFA3                0.157472   0.029590   5.322 1.03e-07 ***
## SMA_RGFA4               -0.180015   0.031830  -5.656 1.55e-08 ***
## SMA_RGFA5                0.005761   0.034115   0.169   0.8659    
## SMA_RGFA6               -0.008236   0.028576  -0.288   0.7732    
## SMA_RGFA7               -0.004237   0.033900  -0.125   0.9005    
## SMA_RGFA8               -0.049918   0.038989  -1.280   0.2004    
## age                      0.004140   0.004919   0.842   0.4000    
## femaleyes               -0.008176   0.077878  -0.105   0.9164    
## race.ethnicityBlack     -0.352280   0.140841  -2.501   0.0124 *  
## race.ethnicityHispanic  -0.507607   0.122478  -4.144 3.41e-05 ***
## race.ethnicityAsian      0.872363   0.183671   4.750 2.04e-06 ***
## race.ethnicityOther     -0.232362   0.113001  -2.056   0.0398 *  
## high.educHS Degree       0.463863   0.255575   1.815   0.0695 .  
## high.educCollege Degree  0.548912   0.266188   2.062   0.0392 *  
## high.educBachelor        0.583931   0.263019   2.220   0.0264 *  
## high.educHigher          0.551252   0.270020   2.042   0.0412 *  
## high.educ_miss           0.448188   0.855908   0.524   0.6005    
## marriedyes              -0.071178   0.088221  -0.807   0.4198    
## HHInc[50K - 100K]       -0.173915   0.103757  -1.676   0.0937 .  
## HHInc[>100K]            -0.074867   0.116141  -0.645   0.5192    
## HHInc_miss              -0.347560   0.165330  -2.102   0.0355 *  
## demo_prnt_age_v2        -0.001006   0.005850  -0.172   0.8634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.34416    0.06106   22.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(kidsi_sma_ZIP,type="est",vline.color = "red",order.terms = c(1:48),show.values = TRUE,value.offset = .75,title = "ZIP Screen Media GFA: Total Youth Suicide Items, standard zero inflation")

# With parameterized zero inflation coefficients:

kidsi_cov_pZIP <- myglmmTMBall(mydepvar,"null",covars,nestvars,glmmadmbdata,"poisson",2)
summary(kidsi_cov_pZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ age + female + race.ethnicity + high.educ + married +  
##     HHInc + demo_prnt_age_v2 + (1 | site_name)
## Zero inflation:              
## ~age + female + race.ethnicity + high.educ + married + HHInc + 
##     demo_prnt_age_v2 + (1 | site_name)
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5282.4   5512.4  -2605.2   5210.4     4359 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.0268   0.1637  
## Number of obs: 4395, groups:  site_name, 20
## 
## Zero-inflation model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.05276  0.2297  
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.241435   0.735876  -1.687  0.09160 .  
## age                      0.016245   0.005523   2.941  0.00327 ** 
## femaleyes               -0.049145   0.079870  -0.615  0.53835    
## race.ethnicityBlack     -0.019083   0.142792  -0.134  0.89369    
## race.ethnicityHispanic  -0.591926   0.143425  -4.127 3.67e-05 ***
## race.ethnicityAsian      0.545136   0.181118   3.010  0.00261 ** 
## race.ethnicityOther     -0.207016   0.129575  -1.598  0.11012    
## high.educHS Degree       0.226943   0.314479   0.722  0.47051    
## high.educCollege Degree  0.428824   0.323927   1.324  0.18556    
## high.educBachelor        0.308806   0.320961   0.962  0.33599    
## high.educHigher          0.120954   0.329141   0.367  0.71326    
## high.educ_miss           0.190259   1.028413   0.185  0.85323    
## marriedyes              -0.182374   0.094151  -1.937  0.05274 .  
## HHInc[50K - 100K]       -0.271821   0.111490  -2.438  0.01477 *  
## HHInc[>100K]            -0.279924   0.121272  -2.308  0.02099 *  
## HHInc_miss              -0.426012   0.194408  -2.191  0.02843 *  
## demo_prnt_age_v2        -0.001207   0.006349  -0.190  0.84923    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.773336   0.963299  -1.841 0.065636 .  
## age                      0.024016   0.007231   3.321 0.000896 ***
## femaleyes                0.418310   0.105092   3.980 6.88e-05 ***
## race.ethnicityBlack      0.294053   0.195439   1.505 0.132432    
## race.ethnicityHispanic  -0.026339   0.181171  -0.145 0.884409    
## race.ethnicityAsian     -0.043778   0.300127  -0.146 0.884027    
## race.ethnicityOther     -0.316557   0.175852  -1.800 0.071839 .  
## high.educHS Degree      -0.068647   0.361937  -0.190 0.849571    
## high.educCollege Degree -0.066718   0.376836  -0.177 0.859470    
## high.educBachelor       -0.109332   0.373002  -0.293 0.769434    
## high.educHigher         -0.194319   0.381802  -0.509 0.610785    
## high.educ_miss          -1.149954   1.596462  -0.720 0.471332    
## marriedyes               0.109657   0.128883   0.851 0.394869    
## HHInc[50K - 100K]       -0.176440   0.158223  -1.115 0.264793    
## HHInc[>100K]            -0.209483   0.174863  -1.198 0.230923    
## HHInc_miss              -0.139088   0.240464  -0.578 0.562985    
## demo_prnt_age_v2         0.012351   0.008494   1.454 0.145916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(kidsi_cov_ZIP,title = "ZIP Covariates: Youth Suicide Items, parameterized standard zero inflation")

kidsi_sma_pZIP <- myglmmTMBall(mydepvar,indepvars,covars,nestvars,glmmadmbdata,"poisson",2)
summary(kidsi_sma_pZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 +  
##     SMA_RGFA5 + SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female +  
##     race.ethnicity + high.educ + married + HHInc + demo_prnt_age_v2 +  
##     (1 | site_name)
## Zero inflation:              
## ~SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 + SMA_RGFA5 + 
##     SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female + race.ethnicity + 
##     high.educ + married + HHInc + demo_prnt_age_v2 + (1 | site_name)
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5086.9   5419.1  -2491.5   4982.9     4343 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.03801  0.195   
## Number of obs: 4395, groups:  site_name, 20
## 
## Zero-inflation model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.08153  0.2855  
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.796595   0.760427  -2.363 0.018147 *  
## SMA_RGFA1                0.150937   0.033159   4.552 5.32e-06 ***
## SMA_RGFA2                0.161095   0.044219   3.643 0.000269 ***
## SMA_RGFA3                0.145464   0.032222   4.514 6.35e-06 ***
## SMA_RGFA4               -0.144378   0.035735  -4.040 5.34e-05 ***
## SMA_RGFA5                0.005826   0.037755   0.154 0.877360    
## SMA_RGFA6               -0.032119   0.031518  -1.019 0.308168    
## SMA_RGFA7               -0.038871   0.038668  -1.005 0.314777    
## SMA_RGFA8               -0.042115   0.042059  -1.001 0.316660    
## age                      0.014251   0.005584   2.552 0.010713 *  
## femaleyes                0.082558   0.086099   0.959 0.337620    
## race.ethnicityBlack     -0.272001   0.153604  -1.771 0.076595 .  
## race.ethnicityHispanic  -0.534914   0.147058  -3.637 0.000275 ***
## race.ethnicityAsian      0.705959   0.186036   3.795 0.000148 ***
## race.ethnicityOther     -0.322950   0.132112  -2.445 0.014505 *  
## high.educHS Degree       0.497293   0.335952   1.480 0.138807    
## high.educCollege Degree  0.588085   0.343281   1.713 0.086688 .  
## high.educBachelor        0.561256   0.342393   1.639 0.101169    
## high.educHigher          0.486957   0.351634   1.385 0.166102    
## high.educ_miss          -0.036011   1.111725  -0.032 0.974160    
## marriedyes              -0.041464   0.099657  -0.416 0.677358    
## HHInc[50K - 100K]       -0.249227   0.113220  -2.201 0.027717 *  
## HHInc[>100K]            -0.189944   0.125993  -1.508 0.131663    
## HHInc_miss              -0.386429   0.199327  -1.939 0.052541 .  
## demo_prnt_age_v2         0.003188   0.006452   0.494 0.621251    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.673567   1.057786  -1.582  0.11362    
## SMA_RGFA1               -0.519031   0.054048  -9.603  < 2e-16 ***
## SMA_RGFA2               -0.047570   0.062658  -0.759  0.44773    
## SMA_RGFA3               -0.104250   0.052059  -2.003  0.04523 *  
## SMA_RGFA4                0.074173   0.053141   1.396  0.16278    
## SMA_RGFA5               -0.002419   0.054169  -0.045  0.96438    
## SMA_RGFA6               -0.006655   0.051490  -0.129  0.89716    
## SMA_RGFA7               -0.037154   0.055772  -0.666  0.50531    
## SMA_RGFA8               -0.036364   0.061159  -0.595  0.55213    
## age                      0.024438   0.007716   3.167  0.00154 ** 
## femaleyes                0.254070   0.118495   2.144  0.03202 *  
## race.ethnicityBlack      0.196644   0.217887   0.903  0.36679    
## race.ethnicityHispanic  -0.093767   0.195009  -0.481  0.63063    
## race.ethnicityAsian     -0.146806   0.307360  -0.478  0.63291    
## race.ethnicityOther     -0.319937   0.191950  -1.667  0.09556 .  
## high.educHS Degree       0.109608   0.427783   0.256  0.79778    
## high.educCollege Degree  0.042918   0.442033   0.097  0.92265    
## high.educBachelor       -0.054545   0.439962  -0.124  0.90133    
## high.educHigher         -0.147065   0.449414  -0.327  0.74349    
## high.educ_miss          -1.023216   2.199266  -0.465  0.64175    
## marriedyes               0.048824   0.139015   0.351  0.72543    
## HHInc[50K - 100K]       -0.310991   0.168712  -1.843  0.06528 .  
## HHInc[>100K]            -0.371627   0.186472  -1.993  0.04627 *  
## HHInc_miss              -0.208031   0.260568  -0.798  0.42465    
## demo_prnt_age_v2         0.013092   0.008973   1.459  0.14453    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(kidsi_sma_ZIP,type="est",vline.color = "red",order.terms = c(1:48),show.values = TRUE,value.offset = .75,title = "ZIP Screen Media GFA: Total Youth Suicide Items, standard zero inflation")

# Negative Binomial Model
# Unfortunately the coefficient plot does not work with nbinom
# could use a self-made plot from here:
# https://github.com/glmmTMB/glmmTMB/blob/master/misc/salamanders.rmd
# https://www.fromthebottomoftheheap.net/2017/05/04/compare-mgcv-with-glmmTMB/
# without zero inflation coefficients:


kidsi_cov_ZINB <- myglmmTMBall(mydepvar,"null",covars,nestvars,glmmadmbdata,"nbinom2",1)
summary(kidsi_cov_ZINB)

##  Family: nbinom2  ( log )
## Formula:          
## kidsitotal ~ age + female + race.ethnicity + high.educ + married +  
##     HHInc + demo_prnt_age_v2 + (1 | site_name)
## Zero inflation:              ~1
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5037.4   5165.2  -2498.7   4997.4     4375 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.04053  0.2013  
## Number of obs: 4395, groups:  site_name, 20
## 
## Overdispersion parameter for nbinom2 family (): 0.132 
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.363242   0.956762  -0.380 0.704199    
## age                     -0.002915   0.007261  -0.401 0.688075    
## femaleyes               -0.372921   0.106119  -3.514 0.000441 ***
## race.ethnicityBlack     -0.250153   0.193782  -1.291 0.196740    
## race.ethnicityHispanic  -0.557152   0.160487  -3.472 0.000517 ***
## race.ethnicityAsian      0.623805   0.328237   1.900 0.057371 .  
## race.ethnicityOther      0.060307   0.178928   0.337 0.736081    
## high.educHS Degree       0.335598   0.325887   1.030 0.303106    
## high.educCollege Degree  0.497772   0.344569   1.445 0.148564    
## high.educBachelor        0.420569   0.335751   1.253 0.210344    
## high.educHigher          0.345124   0.345630   0.999 0.318019    
## high.educ_miss           0.968624   1.602284   0.605 0.545493    
## marriedyes              -0.297675   0.130243  -2.286 0.022282 *  
## HHInc[50K - 100K]       -0.136530   0.157948  -0.864 0.387371    
## HHInc[>100K]            -0.101480   0.171444  -0.592 0.553908    
## HHInc_miss              -0.364589   0.230175  -1.584 0.113202    
## demo_prnt_age_v2        -0.011077   0.008706  -1.272 0.203283    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   -16.39    2128.55  -0.008    0.994

kidsi_sma_ZINB <- myglmmTMBall(mydepvar,indepvars,covars,nestvars,glmmadmbdata,"nbinom2",1)
summary(kidsi_sma_ZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 +  
##     SMA_RGFA5 + SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female +  
##     race.ethnicity + high.educ + married + HHInc + demo_prnt_age_v2 +  
##     (1 | site_name)
## Zero inflation:              ~1
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5183.7   5356.2  -2564.8   5129.7     4368 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.0649   0.2548  
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.631630   0.664875  -0.950   0.3421    
## SMA_RGFA1                0.270817   0.030413   8.905  < 2e-16 ***
## SMA_RGFA2                0.179537   0.039740   4.518 6.25e-06 ***
## SMA_RGFA3                0.157472   0.029590   5.322 1.03e-07 ***
## SMA_RGFA4               -0.180015   0.031830  -5.656 1.55e-08 ***
## SMA_RGFA5                0.005761   0.034115   0.169   0.8659    
## SMA_RGFA6               -0.008236   0.028576  -0.288   0.7732    
## SMA_RGFA7               -0.004237   0.033900  -0.125   0.9005    
## SMA_RGFA8               -0.049918   0.038989  -1.280   0.2004    
## age                      0.004140   0.004919   0.842   0.4000    
## femaleyes               -0.008176   0.077878  -0.105   0.9164    
## race.ethnicityBlack     -0.352280   0.140841  -2.501   0.0124 *  
## race.ethnicityHispanic  -0.507607   0.122478  -4.144 3.41e-05 ***
## race.ethnicityAsian      0.872363   0.183671   4.750 2.04e-06 ***
## race.ethnicityOther     -0.232362   0.113001  -2.056   0.0398 *  
## high.educHS Degree       0.463863   0.255575   1.815   0.0695 .  
## high.educCollege Degree  0.548912   0.266188   2.062   0.0392 *  
## high.educBachelor        0.583931   0.263019   2.220   0.0264 *  
## high.educHigher          0.551252   0.270020   2.042   0.0412 *  
## high.educ_miss           0.448188   0.855908   0.524   0.6005    
## marriedyes              -0.071178   0.088221  -0.807   0.4198    
## HHInc[50K - 100K]       -0.173915   0.103757  -1.676   0.0937 .  
## HHInc[>100K]            -0.074867   0.116141  -0.645   0.5192    
## HHInc_miss              -0.347560   0.165330  -2.102   0.0355 *  
## demo_prnt_age_v2        -0.001006   0.005850  -0.172   0.8634    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.34416    0.06106   22.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

kidsi_sma_pZINB <- myglmmTMBall(mydepvar,indepvars,covars,nestvars,glmmadmbdata,"poisson",2)
summary(kidsi_sma_pZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 +  
##     SMA_RGFA5 + SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female +  
##     race.ethnicity + high.educ + married + HHInc + demo_prnt_age_v2 +  
##     (1 | site_name)
## Zero inflation:              
## ~SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 + SMA_RGFA5 + 
##     SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female + race.ethnicity + 
##     high.educ + married + HHInc + demo_prnt_age_v2 + (1 | site_name)
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   5086.9   5419.1  -2491.5   4982.9     4343 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.03801  0.195   
## Number of obs: 4395, groups:  site_name, 20
## 
## Zero-inflation model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.08153  0.2855  
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.796595   0.760427  -2.363 0.018147 *  
## SMA_RGFA1                0.150937   0.033159   4.552 5.32e-06 ***
## SMA_RGFA2                0.161095   0.044219   3.643 0.000269 ***
## SMA_RGFA3                0.145464   0.032222   4.514 6.35e-06 ***
## SMA_RGFA4               -0.144378   0.035735  -4.040 5.34e-05 ***
## SMA_RGFA5                0.005826   0.037755   0.154 0.877360    
## SMA_RGFA6               -0.032119   0.031518  -1.019 0.308168    
## SMA_RGFA7               -0.038871   0.038668  -1.005 0.314777    
## SMA_RGFA8               -0.042115   0.042059  -1.001 0.316660    
## age                      0.014251   0.005584   2.552 0.010713 *  
## femaleyes                0.082558   0.086099   0.959 0.337620    
## race.ethnicityBlack     -0.272001   0.153604  -1.771 0.076595 .  
## race.ethnicityHispanic  -0.534914   0.147058  -3.637 0.000275 ***
## race.ethnicityAsian      0.705959   0.186036   3.795 0.000148 ***
## race.ethnicityOther     -0.322950   0.132112  -2.445 0.014505 *  
## high.educHS Degree       0.497293   0.335952   1.480 0.138807    
## high.educCollege Degree  0.588085   0.343281   1.713 0.086688 .  
## high.educBachelor        0.561256   0.342393   1.639 0.101169    
## high.educHigher          0.486957   0.351634   1.385 0.166102    
## high.educ_miss          -0.036011   1.111725  -0.032 0.974160    
## marriedyes              -0.041464   0.099657  -0.416 0.677358    
## HHInc[50K - 100K]       -0.249227   0.113220  -2.201 0.027717 *  
## HHInc[>100K]            -0.189944   0.125993  -1.508 0.131663    
## HHInc_miss              -0.386429   0.199327  -1.939 0.052541 .  
## demo_prnt_age_v2         0.003188   0.006452   0.494 0.621251    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.673567   1.057786  -1.582  0.11362    
## SMA_RGFA1               -0.519031   0.054048  -9.603  < 2e-16 ***
## SMA_RGFA2               -0.047570   0.062658  -0.759  0.44773    
## SMA_RGFA3               -0.104250   0.052059  -2.003  0.04523 *  
## SMA_RGFA4                0.074173   0.053141   1.396  0.16278    
## SMA_RGFA5               -0.002419   0.054169  -0.045  0.96438    
## SMA_RGFA6               -0.006655   0.051490  -0.129  0.89716    
## SMA_RGFA7               -0.037154   0.055772  -0.666  0.50531    
## SMA_RGFA8               -0.036364   0.061159  -0.595  0.55213    
## age                      0.024438   0.007716   3.167  0.00154 ** 
## femaleyes                0.254070   0.118495   2.144  0.03202 *  
## race.ethnicityBlack      0.196644   0.217887   0.903  0.36679    
## race.ethnicityHispanic  -0.093767   0.195009  -0.481  0.63063    
## race.ethnicityAsian     -0.146806   0.307360  -0.478  0.63291    
## race.ethnicityOther     -0.319937   0.191950  -1.667  0.09556 .  
## high.educHS Degree       0.109608   0.427783   0.256  0.79778    
## high.educCollege Degree  0.042918   0.442033   0.097  0.92265    
## high.educBachelor       -0.054545   0.439962  -0.124  0.90133    
## high.educHigher         -0.147065   0.449414  -0.327  0.74349    
## high.educ_miss          -1.023216   2.199266  -0.465  0.64175    
## marriedyes               0.048824   0.139015   0.351  0.72543    
## HHInc[50K - 100K]       -0.310991   0.168712  -1.843  0.06528 .  
## HHInc[>100K]            -0.371627   0.186472  -1.993  0.04627 *  
## HHInc_miss              -0.208031   0.260568  -0.798  0.42465    
## demo_prnt_age_v2         0.013092   0.008973   1.459  0.14453    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Non zero-inflated models
kidsi_sma_nZINB <- myglmmTMBall(mydepvar,indepvars,covars,nestvars,glmmadmbdata,"nbinom2",0)
summary(kidsi_sma_nZINB)

##  Family: nbinom2  ( log )
## Formula:          
## kidsitotal ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 +  
##     SMA_RGFA5 + SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female +  
##     race.ethnicity + high.educ + married + HHInc + demo_prnt_age_v2 +  
##     (1 | site_name)
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   4900.0   5072.5  -2423.0   4846.0     4368 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.08234  0.2869  
## Number of obs: 4395, groups:  site_name, 20
## 
## Overdispersion parameter for nbinom2 family (): 0.17 
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.569366   0.944845  -0.603 0.546773    
## SMA_RGFA1                0.567331   0.054598  10.391  < 2e-16 ***
## SMA_RGFA2                0.217421   0.056221   3.867 0.000110 ***
## SMA_RGFA3                0.178284   0.051807   3.441 0.000579 ***
## SMA_RGFA4               -0.213371   0.052503  -4.064 4.82e-05 ***
## SMA_RGFA5                0.002363   0.049728   0.048 0.962104    
## SMA_RGFA6               -0.008862   0.048990  -0.181 0.856458    
## SMA_RGFA7               -0.025933   0.054319  -0.477 0.633069    
## SMA_RGFA8               -0.050206   0.059887  -0.838 0.401836    
## age                     -0.009226   0.007126  -1.295 0.195398    
## femaleyes               -0.133102   0.110337  -1.206 0.227695    
## race.ethnicityBlack     -0.474283   0.202785  -2.339 0.019343 *  
## race.ethnicityHispanic  -0.333063   0.162275  -2.052 0.040125 *  
## race.ethnicityAsian      0.888068   0.305061   2.911 0.003601 ** 
## race.ethnicityOther     -0.032561   0.175051  -0.186 0.852440    
## high.educHS Degree       0.609914   0.325780   1.872 0.061184 .  
## high.educCollege Degree  0.667092   0.343698   1.941 0.052268 .  
## high.educBachelor        0.736175   0.337593   2.181 0.029209 *  
## high.educHigher          0.759336   0.347056   2.188 0.028674 *  
## high.educ_miss           0.775996   1.509043   0.514 0.607091    
## marriedyes              -0.118224   0.127537  -0.927 0.353936    
## HHInc[50K - 100K]       -0.004128   0.157069  -0.026 0.979035    
## HHInc[>100K]             0.142491   0.170313   0.837 0.402794    
## HHInc_miss              -0.173630   0.226275  -0.767 0.442879    
## demo_prnt_age_v2        -0.009669   0.008501  -1.137 0.255345    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(kidsi_sma_nZINB,type="est",vline.color = "red",order.terms = c(1:24),show.values = TRUE,value.offset = .75,title = "nZINB Screen Media GFA: Total Youth Suicide Items, no zero inflation")

kidsi_sma_nZIP <- myglmmTMBall(mydepvar,indepvars,covars,nestvars,glmmadmbdata,"poisson",0)
summary(kidsi_sma_nZIP)

##  Family: poisson  ( log )
## Formula:          
## kidsitotal ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 +  
##     SMA_RGFA5 + SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female +  
##     race.ethnicity + high.educ + married + HHInc + demo_prnt_age_v2 +  
##     (1 | site_name)
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   6265.9   6432.0  -3107.0   6213.9     4369 
## 
## Random effects:
## 
## Conditional model:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.1156   0.34    
## Number of obs: 4395, groups:  site_name, 20
## 
## Conditional model:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -1.300797   0.556971  -2.335  0.01952 *  
## SMA_RGFA1                0.500736   0.024646  20.317  < 2e-16 ***
## SMA_RGFA2                0.182341   0.033653   5.418 6.02e-08 ***
## SMA_RGFA3                0.209126   0.025445   8.219  < 2e-16 ***
## SMA_RGFA4               -0.174227   0.026173  -6.657 2.80e-11 ***
## SMA_RGFA5               -0.014301   0.029327  -0.488  0.62582    
## SMA_RGFA6               -0.007878   0.026359  -0.299  0.76504    
## SMA_RGFA7                0.014081   0.027977   0.503  0.61475    
## SMA_RGFA8               -0.024542   0.032412  -0.757  0.44894    
## age                     -0.001893   0.004067  -0.466  0.64157    
## femaleyes               -0.120460   0.064354  -1.872  0.06123 .  
## race.ethnicityBlack     -0.355091   0.116719  -3.042  0.00235 ** 
## race.ethnicityHispanic  -0.510717   0.105045  -4.862 1.16e-06 ***
## race.ethnicityAsian      0.813728   0.158090   5.147 2.64e-07 ***
## race.ethnicityOther     -0.087106   0.099255  -0.878  0.38016    
## high.educHS Degree       0.298566   0.210386   1.419  0.15586    
## high.educCollege Degree  0.497560   0.218590   2.276  0.02283 *  
## high.educBachelor        0.584017   0.216650   2.696  0.00702 ** 
## high.educHigher          0.533431   0.222081   2.402  0.01631 *  
## high.educ_miss           0.873833   0.749378   1.166  0.24358    
## marriedyes              -0.126554   0.072868  -1.737  0.08243 .  
## HHInc[50K - 100K]       -0.037590   0.089456  -0.420  0.67433    
## HHInc[>100K]             0.036737   0.099917   0.368  0.71312    
## HHInc_miss              -0.309998   0.137875  -2.248  0.02455 *  
## demo_prnt_age_v2        -0.006001   0.004834  -1.241  0.21444    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot_model(kidsi_sma_nZIP,type="est",vline.color = "red",order.terms = c(1:24),show.values = TRUE,value.offset = .75,title = "nZIP Screen Media GFA: Total Youth Suicide Items, no zero inflation")

# Compare the AIC

bictab <- BICtab(kidsi_cov_ZIP,
                 kidsi_sma_ZIP,
                 kidsi_cov_pZIP,
                 kidsi_sma_pZIP,
                 kidsi_cov_ZINB,
                 kidsi_sma_ZINB,
                 kidsi_sma_pZINB,
                 kidsi_sma_nZIP,
                 kidsi_sma_nZINB,logLik=TRUE,base=TRUE)

tabAsStringMatrix <- print(bictab, printToggle = FALSE, noSpaces = TRUE)

##                 logLik  BIC     dLogLik dBIC    df
## kidsi_sma_nZINB -2423.0  5072.5   684.0     0.0 27
## kidsi_sma_ZINB  -2423.0  5080.9   684.0     8.4 28
## kidsi_cov_ZINB  -2498.7  5165.2   608.3    92.7 20
## kidsi_sma_ZIP   -2564.8  5356.2   542.1   283.7 27
## kidsi_cov_ZIP   -2628.9  5417.1   478.1   344.7 19
## kidsi_sma_pZIP  -2491.5  5419.1   615.5   346.6 52
## kidsi_sma_pZINB -2491.5  5419.1   615.5   346.6 52
## kidsi_cov_pZIP  -2605.2  5512.4   501.8   439.9 36
## kidsi_sma_nZIP  -3107.0  6432.0     0.0  1359.6 26

kable(tabAsStringMatrix, "html",caption = "Bayesian Information Criterion: Model Comparison") %>% kable_styling(bootstrap_options = c( "hover", "condensed", "responsive"),font_size = 11)

Bayesian Information Criterion: Model Comparison

	logLik	BIC	dLogLik	dBIC	df
kidsi_sma_nZINB	-2423.0	5072.5	684.0	0.0	27
kidsi_sma_ZINB	-2423.0	5080.9	684.0	8.4	28
kidsi_cov_ZINB	-2498.7	5165.2	608.3	92.7	20
kidsi_sma_ZIP	-2564.8	5356.2	542.1	283.7	27
kidsi_cov_ZIP	-2628.9	5417.1	478.1	344.7	19
kidsi_sma_pZIP	-2491.5	5419.1	615.5	346.6	52
kidsi_sma_pZINB	-2491.5	5419.1	615.5	346.6	52
kidsi_cov_pZIP	-2605.2	5512.4	501.8	439.9	36
kidsi_sma_nZIP	-3107.0	6432.0	0.0	1359.6	26

# Plotting Results

# Remove high variance variables:
myremove <- c(1,20)

# Forming different data frames for each model:
dfmymodel <- createCoeftab(kidsi_sma_ZINB,"ZINB, standard zero inflation",c(GFAselectLabel,colabels),myremove)
dfmymodel2 <- createCoeftab(kidsi_sma_pZINB,"ZINB, parameterized zero inflation",c(GFAselectLabel,colabels),myremove)
dfmymodel3  <- createCoeftab(kidsi_sma_ZIP,"ZIP, standard zero inflation",c(GFAselectLabel,colabels),myremove)
dfmymodel4  <- createCoeftab(kidsi_sma_nZINB,"ZINB, no zero inflation",c(GFAselectLabel,colabels),myremove)

allmodels <- rbind(dfmymodel,dfmymodel2,dfmymodel3,dfmymodel4)

# Plotting all models
ggplot(allmodels, aes(x = estimate, y = term, colour = model, shape = model, xmax = upper, xmin = lower)) + theme_minimal() +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  
    geom_pointrangeh(position = position_dodgev(height = 0.5)) +
    labs(y = NULL,
         x = "Regression estimate +/- CI",
         title = "Zero Inflated Models")

# Exponentiated Model
ggplot(allmodels, aes(x = expestimate, y = term, colour = model, shape = model, xmax = expupper, xmin = explower)) + theme_minimal() +
  geom_vline(xintercept = 1, colour = "grey60", linetype = 2) +
  
    geom_pointrangeh(position = position_dodgev(height = 0.5)) +
    labs(y = NULL,
         x = "Risk Ratio estimates +/- CI",
         title = "Zero Inflated Models")

# Plot best model
ggplot(dfmymodel4, aes(x = estimate, y = term, xmax = upper, xmin = lower)) + theme_minimal() +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  geom_text(aes(label = format(estimate, digits=2, nsmall=2)),nudge_y = -0.3,nudge_x = 0.3,size = 3) +
    geom_pointrangeh(position = position_dodgev(height = 0.5)) +
    labs(y = NULL,
         x = "Regression estimate +/- CI",
         title = "Best Non-Zero Inflated Model: Negative Binomial")

# Interpretation: https://stats.idre.ucla.edu/stata/output/negative-binomial-regression/
# http://fmwww.bc.edu/ec-c/s2013/327/S5CountCategorical0511.slides.pdf
# for a one unit change in the predictor variable, the difference in the logs of expected counts of the response variable is expected to change by the respective regression coefficient, given the other predictor variables in the model are held constant

# Exponentiated Model
ggplot(dfmymodel4, aes(x = expestimate, y = term, xmax = expupper, xmin = explower)) + theme_minimal() +
  geom_vline(xintercept = 1, colour = "grey60", linetype = 2) +
  geom_text(aes(label = format(expestimate, digits=2, nsmall=2)),nudge_y = -0.3,nudge_x = 0.3,size = 3) +
    geom_pointrangeh(position = position_dodgev(height = 0.5)) +
    labs(y = NULL,
         x = "Risk Ratio estimates +/- CI",
         title = "Best Non-Zero Inflated Model: Negative Binomial")

# Interpretation: http://www.mathematica-journal.com/2013/06/negative-binomial-regression/
# Incidence Rate Ratio (IRR) for each variable, which is obtained by exponentiating each coefficient

# Plot a grid of effects:
theme_set(theme_sjplot())
plotlist <- list()
depvars <- paste0("SMA_RGFA",rep(1:4))

# Arrange the marginal plots in a grid of 4 without scatter plots:
for(i in 1:length(depvars)){
  p1 <- plot_model(kidsi_sma_nZINB,type = "pred",terms = c(depvars[i],"female"),title = "",axis.title = c(GFAselectLabel[i],"Youth Total SI Items"))
  plotlist[[i]] <- p1
  }

## Following variables had many unique values and were prettified: SMA_RGFA1. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

## Following variables had many unique values and were prettified: SMA_RGFA2. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

## Following variables had many unique values and were prettified: SMA_RGFA3. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

## Following variables had many unique values and were prettified: SMA_RGFA4. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

mygridtitle <- paste("Marginal Plots: ",sep="")
grid.arrange(plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol = 2,top=mygridtitle)

# Arrange the marginal plots in a grid of 4 with scatter plots:
for(i in 1:length(depvars)){
  p1 <- plot_model(kidsi_sma_nZINB,type = "pred",terms = c(depvars[i],"female"),title = "",axis.title = c(GFAselectLabel[i],"Youth Total SI Items"),show.data = TRUE)
  plotlist[[i]] <- p1
}

## Following variables had many unique values and were prettified: SMA_RGFA1. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

## Following variables had many unique values and were prettified: SMA_RGFA2. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

## Following variables had many unique values and were prettified: SMA_RGFA3. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

## Following variables had many unique values and were prettified: SMA_RGFA4. Use `pretty = FALSE` to get smoother plots with all values, however, at the cost of increased memory usage.

mygridtitle <- paste("Marginal Plots: ",sep="")
grid.arrange(plotlist[[1]],plotlist[[2]],plotlist[[3]],plotlist[[4]],ncol = 2,top=mygridtitle)

This is the logistic regression with both site and family as a random effect. Although there is an output and reasonable stastistics, the routine threw a warning about possible convergence problems.

# Set up the variables:

covars <- c("age","female","race.ethnicity","high.educ","married","HHInc","demo_prnt_age_v2")
colabels <- c("Age","Female","Race: Black","Race: Hispanic","Race: Asian","Race: Other","Parental Education: HS","Parental Education: College","Parental Education: Bachelor","Parental Education: > Bachelor","Parental Education: missing","Married","Household Income: 50-100K","Household Income: > 100K","Household Income: miss","Parental Age")

## Select the nesting variables: site and twin status
nestvars <- c("site_name","FamilyID")
# Independent variables:
indepvars <- paste0("SMA_RGFA",rep(1:8))
# dependent variables:
depvars <- c("KidsSIyes","ParSIyes")

GFAselectLabel <- paste0("SMA_RGFA",rep(1:8))

# Creating a data frame:
glmmadmbdata <- currdata[complete.cases(currdata[,c(covars,nestvars,depvars,indepvars)]),c(covars,nestvars,depvars,indepvars)]

# Visualizing missing data
aggr(currdata[,c(covars,nestvars,depvars,indepvars)],col = c("blue","orange"),sortVars=TRUE,prop = FALSE, numbers = TRUE, combined = TRUE, cex.lab = 0.5, cex.axis =0.5, cex.numbers =0.5)

## 
##  Variables sorted by number of missings: 
##          Variable Count
##          ParSIyes    70
##         KidsSIyes    53
##  demo_prnt_age_v2    43
##    race.ethnicity     2
##               age     0
##            female     0
##         high.educ     0
##           married     0
##             HHInc     0
##         site_name     0
##          FamilyID     0
##         SMA_RGFA1     0
##         SMA_RGFA2     0
##         SMA_RGFA3     0
##         SMA_RGFA4     0
##         SMA_RGFA5     0
##         SMA_RGFA6     0
##         SMA_RGFA7     0
##         SMA_RGFA8     0

mydepvar <- c("KidsSIyes")

# Compare the model with and without the GFAs:
mygamm4base <- myGAMM4Bin(mydepvar,"null",covars,nestvars,glmmadmbdata)
mygamm4GFA <- myGAMM4Bin(mydepvar,indepvars,covars,nestvars,glmmadmbdata)

# Summary of logistic regression
summary(mygamm4base$mer)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## 
##      AIC      BIC   logLik deviance df.resid 
##   2783.7   2905.1  -1372.9   2745.7     4376 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8201 -0.0125 -0.0043 -0.0009  5.8822 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  FamilyID:site_name (Intercept) 1.724e+02 13.130199
##  site_name          (Intercept) 2.125e-06  0.001458
## Number of obs: 4395, groups:  FamilyID:site_name, 3816; site_name, 20
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## X(Intercept)             -29.26396    5.97755  -4.896 9.80e-07 ***
## Xage                       0.21379    0.03727   5.737 9.66e-09 ***
## Xfemaleyes                -1.40832    0.65128  -2.162 0.030587 *  
## Xrace.ethnicityBlack      -8.87263    2.28436  -3.884 0.000103 ***
## Xrace.ethnicityHispanic   -2.49010    1.44891  -1.719 0.085688 .  
## Xrace.ethnicityAsian       0.17799    3.28808   0.054 0.956830    
## Xrace.ethnicityOther       2.55967    1.29133   1.982 0.047456 *  
## Xhigh.educHS Degree       -2.19768    2.59444  -0.847 0.396953    
## Xhigh.educCollege Degree  -1.18882    2.81928  -0.422 0.673261    
## Xhigh.educBachelor        -0.86491    2.74201  -0.315 0.752435    
## Xhigh.educHigher          -0.66527    2.80094  -0.238 0.812257    
## Xhigh.educ_miss            4.26981    9.68121   0.441 0.659183    
## Xmarriedyes               -0.39745    1.24628  -0.319 0.749794    
## XHHInc[50K - 100K]        -3.00808    1.49577  -2.011 0.044320 *  
## XHHInc[>100K]             -3.32062    1.58725  -2.092 0.036433 *  
## XHHInc_miss               -3.96844    2.33582  -1.699 0.089328 .  
## Xdemo_prnt_age_v2         -0.03457    0.08096  -0.427 0.669324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it

## convergence code: 0
## failure to converge in 10000 evaluations

summary(mygamm4base$gam)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## KidsSIyes ~ age + female + race.ethnicity + high.educ + married + 
##     HHInc + demo_prnt_age_v2
## <environment: 0x7fd543f3bc38>
## 
## Parametric coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -29.26396    5.97755  -4.896 9.80e-07 ***
## age                       0.21379    0.03727   5.737 9.66e-09 ***
## femaleyes                -1.40832    0.65128  -2.162 0.030587 *  
## race.ethnicityBlack      -8.87263    2.28436  -3.884 0.000103 ***
## race.ethnicityHispanic   -2.49010    1.44891  -1.719 0.085688 .  
## race.ethnicityAsian       0.17799    3.28808   0.054 0.956830    
## race.ethnicityOther       2.55967    1.29133   1.982 0.047456 *  
## high.educHS Degree       -2.19768    2.59444  -0.847 0.396953    
## high.educCollege Degree  -1.18882    2.81928  -0.422 0.673261    
## high.educBachelor        -0.86491    2.74201  -0.315 0.752435    
## high.educHigher          -0.66527    2.80094  -0.238 0.812257    
## high.educ_miss            4.26981    9.68121   0.441 0.659183    
## marriedyes               -0.39745    1.24628  -0.319 0.749794    
## HHInc[50K - 100K]        -3.00808    1.49577  -2.011 0.044320 *  
## HHInc[>100K]             -3.32062    1.58725  -2.092 0.036433 *  
## HHInc_miss               -3.96844    2.33582  -1.699 0.089328 .  
## demo_prnt_age_v2         -0.03457    0.08096  -0.427 0.669324    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  -0.00425   
## glmer.ML = 419.44  Scale est. = 1         n = 4395

summary(mygamm4GFA$mer)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## 
##      AIC      BIC   logLik deviance df.resid 
##   3268.2   3440.7  -1607.1   3214.2     4368 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4291 -0.3682 -0.2990 -0.2355  4.7467 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  FamilyID:site_name (Intercept) 0.58233  0.7631  
##  site_name          (Intercept) 0.05495  0.2344  
## Number of obs: 4395, groups:  FamilyID:site_name, 3816; site_name, 20
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## X(Intercept)             -0.230776   1.017997  -0.227 0.820660    
## XSMA_RGFA1                0.590580   0.088840   6.648 2.98e-11 ***
## XSMA_RGFA2                0.121004   0.059190   2.044 0.040922 *  
## XSMA_RGFA3                0.170115   0.051041   3.333 0.000859 ***
## XSMA_RGFA4               -0.148340   0.053008  -2.798 0.005135 ** 
## XSMA_RGFA5                0.002801   0.051022   0.055 0.956218    
## XSMA_RGFA6               -0.017267   0.049750  -0.347 0.728529    
## XSMA_RGFA7                0.018334   0.051260   0.358 0.720595    
## XSMA_RGFA8                0.013130   0.057281   0.229 0.818697    
## Xage                     -0.014486   0.006919  -2.094 0.036287 *  
## Xfemaleyes               -0.226142   0.111591  -2.027 0.042711 *  
## Xrace.ethnicityBlack     -0.371209   0.251975  -1.473 0.140698    
## Xrace.ethnicityHispanic  -0.250511   0.172323  -1.454 0.146022    
## Xrace.ethnicityAsian      0.415554   0.322224   1.290 0.197175    
## Xrace.ethnicityOther      0.143342   0.175433   0.817 0.413883    
## Xhigh.educHS Degree       0.118232   0.317668   0.372 0.709753    
## Xhigh.educCollege Degree  0.222132   0.334909   0.663 0.507164    
## Xhigh.educBachelor        0.285381   0.333071   0.857 0.391547    
## Xhigh.educHigher          0.331531   0.344813   0.961 0.336310    
## Xhigh.educ_miss           0.415166   1.401604   0.296 0.767072    
## Xmarriedyes              -0.068178   0.129646  -0.526 0.598975    
## XHHInc[50K - 100K]        0.174714   0.158230   1.104 0.269516    
## XHHInc[>100K]             0.278856   0.185834   1.501 0.133469    
## XHHInc_miss              -0.014378   0.222728  -0.065 0.948530    
## Xdemo_prnt_age_v2        -0.011587   0.008313  -1.394 0.163399    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 25 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it

## convergence code: 0
## failure to converge in 10000 evaluations

summary(mygamm4GFA$gam)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## KidsSIyes ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 + SMA_RGFA5 + 
##     SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female + race.ethnicity + 
##     high.educ + married + HHInc + demo_prnt_age_v2
## <environment: 0x7fd529ac4ff0>
## 
## Parametric coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.230776   0.915068  -0.252 0.800890    
## SMA_RGFA1                0.590580   0.048544  12.166  < 2e-16 ***
## SMA_RGFA2                0.121004   0.056579   2.139 0.032462 *  
## SMA_RGFA3                0.170115   0.048760   3.489 0.000485 ***
## SMA_RGFA4               -0.148340   0.049693  -2.985 0.002835 ** 
## SMA_RGFA5                0.002801   0.049283   0.057 0.954674    
## SMA_RGFA6               -0.017267   0.047178  -0.366 0.714363    
## SMA_RGFA7                0.018334   0.050364   0.364 0.715836    
## SMA_RGFA8                0.013130   0.056663   0.232 0.816754    
## age                     -0.014486   0.006885  -2.104 0.035395 *  
## femaleyes               -0.226142   0.107885  -2.096 0.036069 *  
## race.ethnicityBlack     -0.371209   0.204715  -1.813 0.069786 .  
## race.ethnicityHispanic  -0.250511   0.158790  -1.578 0.114654    
## race.ethnicityAsian      0.415554   0.315391   1.318 0.187643    
## race.ethnicityOther      0.143342   0.165558   0.866 0.386593    
## high.educHS Degree       0.118232   0.309955   0.381 0.702869    
## high.educCollege Degree  0.222132   0.328777   0.676 0.499275    
## high.educBachelor        0.285381   0.323032   0.883 0.376997    
## high.educHigher          0.331531   0.331524   1.000 0.317300    
## high.educ_miss           0.415166   1.361234   0.305 0.760372    
## marriedyes              -0.068178   0.126371  -0.540 0.589539    
## HHInc[50K - 100K]        0.174714   0.157043   1.113 0.265913    
## HHInc[>100K]             0.278856   0.173978   1.603 0.108974    
## HHInc_miss              -0.014378   0.222246  -0.065 0.948419    
## demo_prnt_age_v2        -0.011587   0.008254  -1.404 0.160374    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.0453   
## glmer.ML = 2698.3  Scale est. = 1         n = 4395

# Visualize the GAMM4 Coefficients:
# https://cran.r-project.org/web/packages/merTools/vignettes/merToolsIntro.html

feEx <- FEsim(mygamm4GFA$mer,1000)
cbind(feEx[,1],round(feEx[,2:4],3))

##                   feEx[, 1]   mean median    sd
## 1              X(Intercept) -0.256 -0.295 0.902
## 2                XSMA_RGFA1  0.593  0.592 0.051
## 3                XSMA_RGFA2  0.121  0.120 0.057
## 4                XSMA_RGFA3  0.170  0.169 0.049
## 5                XSMA_RGFA4 -0.149 -0.148 0.048
## 6                XSMA_RGFA5  0.003  0.002 0.049
## 7                XSMA_RGFA6 -0.021 -0.023 0.048
## 8                XSMA_RGFA7  0.019  0.020 0.049
## 9                XSMA_RGFA8  0.012  0.011 0.058
## 10                     Xage -0.015 -0.015 0.007
## 11               Xfemaleyes -0.225 -0.226 0.109
## 12     Xrace.ethnicityBlack -0.371 -0.374 0.209
## 13  Xrace.ethnicityHispanic -0.250 -0.253 0.156
## 14     Xrace.ethnicityAsian  0.425  0.428 0.306
## 15     Xrace.ethnicityOther  0.141  0.139 0.166
## 16      Xhigh.educHS Degree  0.132  0.134 0.317
## 17 Xhigh.educCollege Degree  0.250  0.260 0.327
## 18       Xhigh.educBachelor  0.309  0.312 0.329
## 19         Xhigh.educHigher  0.352  0.367 0.344
## 20          Xhigh.educ_miss  0.440  0.490 1.376
## 21              Xmarriedyes -0.068 -0.069 0.127
## 22       XHHInc[50K - 100K]  0.171  0.175 0.156
## 23            XHHInc[>100K]  0.274  0.275 0.178
## 24              XHHInc_miss -0.020 -0.013 0.216
## 25        Xdemo_prnt_age_v2 -0.011 -0.012 0.008

feEx$term <- c("Intercept",indepvars,colabels)
feEx$term <- factor(feEx$term,levels = c("Intercept",indepvars,colabels))

# reeduce the data for plotting:
reddata <- feEx[feEx$term!= "Intercept" & feEx$term!= "Parental Education: missing", ]
reddata <- droplevels(reddata)

# theme_bw() + 
gamm4coeff <- ggplot(reddata) + 
  theme_minimal() +
  aes(x = term, ymin = median - 1.96 * sd, 
      ymax = median + 1.96 * sd, y = median) + 
  geom_pointrange() + 
  scale_x_discrete(limits = rev(levels(reddata$term))) +
  geom_text(aes(label = sprintf("%0.2f", round(median, digits = 2))),
                position=position_dodge(width=0.9), vjust=-0.75) +
  geom_hline(yintercept = 0, size = I(1.1), color = I("red")) + 
  coord_flip() + 
  labs(title = paste("Kids SI Logistic Regression",": Median Effect Size",sep=""), 
                    x = "Variables", y = "Standardized Coefficients")

print(gamm4coeff)

ExpfeEx <- data.frame(exp(feEx$mean),exp(feEx$median),exp(feEx$median-1.96*feEx$sd),exp(feEx$median+1.96*feEx$sd))

ExpfeEx <- data.frame(cbind(feEx$term,ExpfeEx))

colnames(ExpfeEx) <- c("term","mean","median","Lower_CI","Upper_CI")
ExpfeEx$term = factor(ExpfeEx$term,levels = c("Intercept",indepvars,colabels))

reddata <- ExpfeEx[ExpfeEx$term!= "Intercept" & ExpfeEx$term!= "Parental Education: missing", ]
reddata <- droplevels(reddata)

# theme_bw() +                      
# Exponentiated Results:
gamm4coeff <- ggplot(reddata) + 
  theme_minimal() +
  aes(x = term, ymin = Lower_CI, 
      ymax = Upper_CI, y = median) + 
  geom_pointrange() + 
  scale_x_discrete(limits = rev(levels(reddata$term))) +
  geom_text(aes(label = sprintf("%0.2f", round(median, digits = 2))),
                position=position_dodge(width=0.9), vjust=-0.75) +
  geom_hline(yintercept = 1, size = I(1.1), color = I("red")) + 
  coord_flip() + 
   labs(title = paste("Kids Sucidal Ideation",": Odds Ratios",sep=""), 
                    x = "Variables", y = "Risk Ratio estimates +/- CI")

print(gamm4coeff)

This is the logistic Regression with site only as the random effect because there were some convergence problems when including both site and family as a random effect.

# Set up the variables:

covars <- c("age","female","race.ethnicity","high.educ","married","HHInc","demo_prnt_age_v2")
colabels <- c("Age","Female","Race: Black","Race: Hispanic","Race: Asian","Race: Other","Parental Education: HS","Parental Education: College","Parental Education: Bachelor","Parental Education: > Bachelor","Parental Education: missing","Married","Household Income: 50-100K","Household Income: > 100K","Household Income: miss","Parental Age")

## Select the nesting variables: site and twin status
nestvars <- c("site_name")
# Independent variables:
indepvars <- paste0("SMA_RGFA",rep(1:8))
# dependent variables:
depvars <- c("KidsSIyes","ParSIyes")

GFAselectLabel <- paste0("SMA_RGFA",rep(1:8))

# Creating a data frame:
glmmadmbdata <- currdata[complete.cases(currdata[,c(covars,nestvars,depvars,indepvars)]),c(covars,nestvars,depvars,indepvars)]

# Visualizing missing data
aggr(currdata[,c(covars,nestvars,depvars,indepvars)],col = c("blue","orange"),sortVars=TRUE,prop = FALSE, numbers = TRUE, combined = TRUE, cex.lab = 0.5, cex.axis =0.5, cex.numbers =0.5)

## 
##  Variables sorted by number of missings: 
##          Variable Count
##          ParSIyes    70
##         KidsSIyes    53
##  demo_prnt_age_v2    43
##    race.ethnicity     2
##               age     0
##            female     0
##         high.educ     0
##           married     0
##             HHInc     0
##         site_name     0
##         SMA_RGFA1     0
##         SMA_RGFA2     0
##         SMA_RGFA3     0
##         SMA_RGFA4     0
##         SMA_RGFA5     0
##         SMA_RGFA6     0
##         SMA_RGFA7     0
##         SMA_RGFA8     0

mydepvar <- c("KidsSIyes")

# Compare the model with and without the GFAs:
mygamm4base <- myGAMM4Bin(mydepvar,"null",covars,nestvars,glmmadmbdata)
mygamm4GFA <- myGAMM4Bin(mydepvar,indepvars,covars,nestvars,glmmadmbdata)

# Summary of logistic regression
summary(mygamm4base$mer)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## 
##      AIC      BIC   logLik deviance df.resid 
##   3420.8   3535.8  -1692.4   3384.8     4377 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.6661 -0.4180 -0.3647 -0.3125  3.9438 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.05065  0.225   
## Number of obs: 4395, groups:  site_name, 20
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## X(Intercept)              0.469804   0.822880   0.571   0.5680    
## Xage                     -0.014858   0.006229  -2.385   0.0171 *  
## Xfemaleyes               -0.427060   0.091935  -4.645  3.4e-06 ***
## Xrace.ethnicityBlack     -0.306963   0.180590  -1.700   0.0892 .  
## Xrace.ethnicityHispanic  -0.279086   0.145618  -1.917   0.0553 .  
## Xrace.ethnicityAsian      0.257080   0.286187   0.898   0.3690    
## Xrace.ethnicityOther      0.194162   0.148827   1.305   0.1920    
## Xhigh.educHS Degree       0.178690   0.285970   0.625   0.5321    
## Xhigh.educCollege Degree  0.266027   0.302290   0.880   0.3788    
## Xhigh.educBachelor        0.249574   0.296579   0.842   0.4001    
## Xhigh.educHigher          0.234990   0.302710   0.776   0.4376    
## Xhigh.educ_miss           1.087923   1.202850   0.904   0.3658    
## Xmarriedyes              -0.187883   0.113537  -1.655   0.0980 .  
## XHHInc[50K - 100K]        0.048831   0.141481   0.345   0.7300    
## XHHInc[>100K]             0.073596   0.156544   0.470   0.6383    
## XHHInc_miss              -0.062693   0.199541  -0.314   0.7534    
## Xdemo_prnt_age_v2        -0.012234   0.007516  -1.628   0.1036    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it

summary(mygamm4base$gam)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## KidsSIyes ~ age + female + race.ethnicity + high.educ + married + 
##     HHInc + demo_prnt_age_v2
## <environment: 0x7fd53e7147c8>
## 
## Parametric coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              0.469804   0.824131   0.570   0.5686    
## age                     -0.014858   0.006237  -2.382   0.0172 *  
## femaleyes               -0.427060   0.092044  -4.640 3.49e-06 ***
## race.ethnicityBlack     -0.306963   0.180731  -1.698   0.0894 .  
## race.ethnicityHispanic  -0.279086   0.145085  -1.924   0.0544 .  
## race.ethnicityAsian      0.257080   0.286354   0.898   0.3693    
## race.ethnicityOther      0.194162   0.148971   1.303   0.1925    
## high.educHS Degree       0.178690   0.286485   0.624   0.5328    
## high.educCollege Degree  0.266027   0.302858   0.878   0.3797    
## high.educBachelor        0.249574   0.297047   0.840   0.4008    
## high.educHigher          0.234990   0.303001   0.776   0.4380    
## high.educ_miss           1.087923   1.202303   0.905   0.3655    
## marriedyes              -0.187883   0.113661  -1.653   0.0983 .  
## HHInc[50K - 100K]        0.048831   0.141698   0.345   0.7304    
## HHInc[>100K]             0.073596   0.156758   0.469   0.6387    
## HHInc_miss              -0.062693   0.199888  -0.314   0.7538    
## demo_prnt_age_v2        -0.012234   0.007528  -1.625   0.1042    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.00727   
## glmer.ML = 3358.5  Scale est. = 1         n = 4395

summary(mygamm4GFA$mer)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## 
##      AIC      BIC   logLik deviance df.resid 
##   3270.5   3436.6  -1609.3   3218.5     4369 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5863 -0.4034 -0.3267 -0.2551  5.6666 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  site_name (Intercept) 0.08448  0.2907  
## Number of obs: 4395, groups:  site_name, 20
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## X(Intercept)             -0.105069   0.859948  -0.122 0.902756    
## XSMA_RGFA1                0.551162   0.045478  12.119  < 2e-16 ***
## XSMA_RGFA2                0.126897   0.052923   2.398 0.016495 *  
## XSMA_RGFA3                0.173123   0.045319   3.820 0.000133 ***
## XSMA_RGFA4               -0.139466   0.046120  -3.024 0.002495 ** 
## XSMA_RGFA5                0.006997   0.046122   0.152 0.879416    
## XSMA_RGFA6               -0.009551   0.044447  -0.215 0.829864    
## XSMA_RGFA7                0.019842   0.047163   0.421 0.673969    
## XSMA_RGFA8                0.015024   0.052976   0.284 0.776725    
## Xage                     -0.015173   0.006445  -2.354 0.018568 *  
## Xfemaleyes               -0.208093   0.100198  -2.077 0.037818 *  
## Xrace.ethnicityBlack     -0.340183   0.190411  -1.787 0.074007 .  
## Xrace.ethnicityHispanic  -0.218884   0.150043  -1.459 0.144618    
## Xrace.ethnicityAsian      0.432521   0.290193   1.490 0.136104    
## Xrace.ethnicityOther      0.114887   0.154229   0.745 0.456325    
## Xhigh.educHS Degree       0.212306   0.296457   0.716 0.473903    
## Xhigh.educCollege Degree  0.320058   0.313358   1.021 0.307074    
## Xhigh.educBachelor        0.388491   0.308335   1.260 0.207682    
## Xhigh.educHigher          0.438993   0.316224   1.388 0.165066    
## Xhigh.educ_miss           0.711396   1.215506   0.585 0.558368    
## Xmarriedyes              -0.057951   0.117264  -0.494 0.621168    
## XHHInc[50K - 100K]        0.167292   0.146085   1.145 0.252137    
## XHHInc[>100K]             0.251290   0.162283   1.548 0.121510    
## XHHInc_miss               0.004942   0.206253   0.024 0.980885    
## Xdemo_prnt_age_v2        -0.010365   0.007687  -1.348 0.177543    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 25 > 12.
## Use print(x, correlation=TRUE)  or
##   vcov(x)     if you need it

summary(mygamm4GFA$gam)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## KidsSIyes ~ SMA_RGFA1 + SMA_RGFA2 + SMA_RGFA3 + SMA_RGFA4 + SMA_RGFA5 + 
##     SMA_RGFA6 + SMA_RGFA7 + SMA_RGFA8 + age + female + race.ethnicity + 
##     high.educ + married + HHInc + demo_prnt_age_v2
## <environment: 0x7fd5325bec00>
## 
## Parametric coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             -0.105069   0.860465  -0.122 0.902814    
## SMA_RGFA1                0.551162   0.045238  12.184  < 2e-16 ***
## SMA_RGFA2                0.126897   0.053001   2.394 0.016655 *  
## SMA_RGFA3                0.173123   0.045270   3.824 0.000131 ***
## SMA_RGFA4               -0.139466   0.046141  -3.023 0.002506 ** 
## SMA_RGFA5                0.006997   0.046175   0.152 0.879552    
## SMA_RGFA6               -0.009551   0.044460  -0.215 0.829911    
## SMA_RGFA7                0.019842   0.047195   0.420 0.674177    
## SMA_RGFA8                0.015024   0.053037   0.283 0.776972    
## age                     -0.015173   0.006449  -2.353 0.018638 *  
## femaleyes               -0.208093   0.100368  -2.073 0.038144 *  
## race.ethnicityBlack     -0.340183   0.190262  -1.788 0.073780 .  
## race.ethnicityHispanic  -0.218884   0.150034  -1.459 0.144593    
## race.ethnicityAsian      0.432521   0.290535   1.489 0.136565    
## race.ethnicityOther      0.114887   0.154395   0.744 0.456809    
## high.educHS Degree       0.212306   0.296981   0.715 0.474684    
## high.educCollege Degree  0.320058   0.313900   1.020 0.307910    
## high.educBachelor        0.388491   0.308789   1.258 0.208352    
## high.educHigher          0.438993   0.316621   1.386 0.165596    
## high.educ_miss           0.711396   1.214148   0.586 0.557928    
## marriedyes              -0.057951   0.117424  -0.494 0.621644    
## HHInc[50K - 100K]        0.167292   0.146307   1.143 0.252857    
## HHInc[>100K]             0.251290   0.162452   1.547 0.121897    
## HHInc_miss               0.004942   0.206620   0.024 0.980919    
## demo_prnt_age_v2        -0.010365   0.007697  -1.347 0.178111    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## R-sq.(adj) =  0.046   
## glmer.ML = 3184.8  Scale est. = 1         n = 4395

# Visualize the GAMM4 Coefficients:
# https://cran.r-project.org/web/packages/merTools/vignettes/merToolsIntro.html

feEx <- FEsim(mygamm4GFA$mer,1000)
cbind(feEx[,1],round(feEx[,2:4],3))

##                   feEx[, 1]   mean median    sd
## 1              X(Intercept) -0.101 -0.105 0.814
## 2                XSMA_RGFA1  0.555  0.556 0.046
## 3                XSMA_RGFA2  0.128  0.129 0.051
## 4                XSMA_RGFA3  0.173  0.174 0.045
## 5                XSMA_RGFA4 -0.138 -0.137 0.046
## 6                XSMA_RGFA5  0.005  0.006 0.046
## 7                XSMA_RGFA6 -0.009 -0.010 0.046
## 8                XSMA_RGFA7  0.021  0.021 0.046
## 9                XSMA_RGFA8  0.014  0.016 0.055
## 10                     Xage -0.015 -0.015 0.006
## 11               Xfemaleyes -0.215 -0.213 0.098
## 12     Xrace.ethnicityBlack -0.349 -0.350 0.195
## 13  Xrace.ethnicityHispanic -0.226 -0.223 0.147
## 14     Xrace.ethnicityAsian  0.423  0.413 0.283
## 15     Xrace.ethnicityOther  0.118  0.113 0.155
## 16      Xhigh.educHS Degree  0.214  0.211 0.289
## 17 Xhigh.educCollege Degree  0.317  0.312 0.302
## 18       Xhigh.educBachelor  0.391  0.398 0.304
## 19         Xhigh.educHigher  0.441  0.446 0.316
## 20          Xhigh.educ_miss  0.747  0.739 1.239
## 21              Xmarriedyes -0.053 -0.053 0.116
## 22       XHHInc[50K - 100K]  0.159  0.162 0.146
## 23            XHHInc[>100K]  0.239  0.236 0.164
## 24              XHHInc_miss  0.011  0.013 0.203
## 25        Xdemo_prnt_age_v2 -0.010 -0.010 0.007

feEx$term <- c("Intercept",indepvars,colabels)
feEx$term <- factor(feEx$term,levels = c("Intercept",indepvars,colabels))

# reeduce the data for plotting:
reddata <- feEx[feEx$term!= "Intercept" & feEx$term!= "Parental Education: missing", ]
reddata <- droplevels(reddata)

# theme_bw() + 
gamm4coeff <- ggplot(reddata) + 
  theme_minimal() +
  aes(x = term, ymin = median - 1.96 * sd, 
      ymax = median + 1.96 * sd, y = median) + 
  geom_pointrange() + 
  scale_x_discrete(limits = rev(levels(reddata$term))) +
  geom_text(aes(label = sprintf("%0.2f", round(median, digits = 2))),
                position=position_dodge(width=0.9), vjust=-0.75) +
  geom_hline(yintercept = 0, size = I(1.1), color = I("red")) + 
  coord_flip() + 
  labs(title = paste("Kids SI Logistic Regression",": Median Effect Size",sep=""), 
                    x = "Variables", y = "Standardized Coefficients")

print(gamm4coeff)

ExpfeEx <- data.frame(exp(feEx$mean),exp(feEx$median),exp(feEx$median-1.96*feEx$sd),exp(feEx$median+1.96*feEx$sd))

ExpfeEx <- data.frame(cbind(feEx$term,ExpfeEx))

colnames(ExpfeEx) <- c("term","mean","median","Lower_CI","Upper_CI")
ExpfeEx$term = factor(ExpfeEx$term,levels = c("Intercept",indepvars,colabels))

reddata <- ExpfeEx[ExpfeEx$term!= "Intercept" & ExpfeEx$term!= "Parental Education: missing", ]
reddata <- droplevels(reddata)

# theme_bw() +                      
# Exponentiated Results:
gamm4coeff <- ggplot(reddata) + 
  theme_minimal() +
  aes(x = term, ymin = Lower_CI, 
      ymax = Upper_CI, y = median) + 
  geom_pointrange() + 
  scale_x_discrete(limits = rev(levels(reddata$term))) +
  geom_text(aes(label = sprintf("%0.2f", round(median, digits = 2))),
                position=position_dodge(width=0.9), vjust=-0.75) +
  geom_hline(yintercept = 1, size = I(1.1), color = I("red")) + 
  coord_flip() + 
   labs(title = paste("Kids Sucidal Ideation",": Odds Ratios",sep=""), 
                    x = "Variables", y = "Risk Ratio estimates +/- CI")

print(gamm4coeff)

These are bar plots of the SMA variables by GFA standard deviations

# This is based on tabledata, which was generated by the standard deviation variable
# of the GFAs
# Pick the color based on: https://htmlcolorcodes.com/color-picker/

tabledata <- currdata[,c(mysdGFA,listvars)]

tablenames <- names(tabledata)
mysmas <- tablenames[grep("_y_",tablenames)]
mysmalabels <- c("Watching TV/Movies per week [h]","Watching Videos per week [h]","Playing Games per week [h]","Texting per week [h]","Social Networking per week [h]","Chatting per week [h]")

for (i in 1:8){
  p <- list()
  for (j in 1: length(mysmas)){
  myplotGFA <- paste0("sdGFA",i)
completetabledata <- tabledata[complete.cases(tabledata[,c(myplotGFA)]),]

p[[j]] <-ggplot(data=completetabledata, aes_string(x=myplotGFA, y=mysmas[j])) +
  geom_bar(stat = "summary", fun.y = "mean",fill="#FF3396") + coord_flip() + 
  # ggtitle("Media Activity and GFA") +
  ylim(0, 22) +
  stat_summary(aes(label=round(..y..,2)), fun.y=mean, geom="text", size=3,
             color = "black",hjust=-0.1) +
  xlab(paste0("GFA",i)) + ylab(mysmalabels[j])  + theme_minimal()
# print(p[[i]])
  }
  do.call("grid.arrange", c(p,list(ncol=2, top=textGrob(paste0("GFA ",i)))))
# do.call(grid.arrange,c(p,main="text"))
}

## Warning: Removed 397 rows containing non-finite values (stat_summary).

## Warning: Removed 397 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 353 rows containing non-finite values (stat_summary).

## Warning: Removed 353 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 342 rows containing non-finite values (stat_summary).

## Warning: Removed 342 rows containing non-finite values (stat_summary).

## Warning: Removed 358 rows containing non-finite values (stat_summary).

## Warning: Removed 358 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 374 rows containing non-finite values (stat_summary).

## Warning: Removed 374 rows containing non-finite values (stat_summary).

## Warning: Removed 307 rows containing non-finite values (stat_summary).

## Warning: Removed 307 rows containing non-finite values (stat_summary).

## Warning: Removed 327 rows containing non-finite values (stat_summary).

## Warning: Removed 327 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 6 rows containing non-finite values (stat_summary).

## Warning: Removed 6 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 354 rows containing non-finite values (stat_summary).

## Warning: Removed 354 rows containing non-finite values (stat_summary).

## Warning: Removed 292 rows containing non-finite values (stat_summary).

## Warning: Removed 292 rows containing non-finite values (stat_summary).

## Warning: Removed 297 rows containing non-finite values (stat_summary).

## Warning: Removed 297 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 401 rows containing non-finite values (stat_summary).

## Warning: Removed 401 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 355 rows containing non-finite values (stat_summary).

## Warning: Removed 355 rows containing non-finite values (stat_summary).

## Warning: Removed 46 rows containing non-finite values (stat_summary).

## Warning: Removed 46 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 36 rows containing non-finite values (stat_summary).

## Warning: Removed 36 rows containing non-finite values (stat_summary).

## Warning: Removed 391 rows containing non-finite values (stat_summary).

## Warning: Removed 391 rows containing non-finite values (stat_summary).

## Warning: Removed 337 rows containing non-finite values (stat_summary).

## Warning: Removed 337 rows containing non-finite values (stat_summary).

## Warning: Removed 352 rows containing non-finite values (stat_summary).

## Warning: Removed 352 rows containing non-finite values (stat_summary).

## Warning: Removed 45 rows containing non-finite values (stat_summary).

## Warning: Removed 45 rows containing non-finite values (stat_summary).

## Warning: Removed 20 rows containing non-finite values (stat_summary).

## Warning: Removed 20 rows containing non-finite values (stat_summary).

## Warning: Removed 34 rows containing non-finite values (stat_summary).

## Warning: Removed 34 rows containing non-finite values (stat_summary).

## Warning: Removed 399 rows containing non-finite values (stat_summary).

## Warning: Removed 399 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 352 rows containing non-finite values (stat_summary).

## Warning: Removed 352 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 389 rows containing non-finite values (stat_summary).

## Warning: Removed 389 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 356 rows containing non-finite values (stat_summary).

## Warning: Removed 356 rows containing non-finite values (stat_summary).

## Warning: Removed 44 rows containing non-finite values (stat_summary).

## Warning: Removed 44 rows containing non-finite values (stat_summary).

## Warning: Removed 21 rows containing non-finite values (stat_summary).

## Warning: Removed 21 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

for (j in 1: length(mysmas)){
p <- list()
for (i in 1:4){
  myplotGFA <- paste0("sdGFA",i)
completetabledata <- tabledata[complete.cases(tabledata[,c(myplotGFA)]),]

p[[i]] <-ggplot(data=completetabledata, aes_string(x=myplotGFA, y=mysmas[j])) +
  geom_bar(stat = "summary", fun.y = "mean",fill="skyblue") + coord_flip() + 
  # ggtitle("Media Activity and GFA") +
  ylim(0, 22) +
  stat_summary(aes(label=round(..y..,2)), fun.y=mean, geom="text", size=3,
             color = "black",hjust=-0.1) +
  xlab(paste0("GFA",i)) + ylab(mysmalabels[j])  + theme_minimal()
# print(p[[i]])
}
do.call("grid.arrange", c(p,list(ncol=2, top=textGrob(paste0("GFA Groups and ",mysmalabels[j])))))
# do.call(grid.arrange,p)
}

## Warning: Removed 397 rows containing non-finite values (stat_summary).

## Warning: Removed 397 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 374 rows containing non-finite values (stat_summary).

## Warning: Removed 374 rows containing non-finite values (stat_summary).

## Warning: Removed 354 rows containing non-finite values (stat_summary).

## Warning: Removed 354 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 339 rows containing non-finite values (stat_summary).

## Warning: Removed 342 rows containing non-finite values (stat_summary).

## Warning: Removed 342 rows containing non-finite values (stat_summary).

## Warning: Removed 307 rows containing non-finite values (stat_summary).

## Warning: Removed 307 rows containing non-finite values (stat_summary).

## Warning: Removed 292 rows containing non-finite values (stat_summary).

## Warning: Removed 292 rows containing non-finite values (stat_summary).

## Warning: Removed 353 rows containing non-finite values (stat_summary).

## Warning: Removed 353 rows containing non-finite values (stat_summary).

## Warning: Removed 358 rows containing non-finite values (stat_summary).

## Warning: Removed 358 rows containing non-finite values (stat_summary).

## Warning: Removed 327 rows containing non-finite values (stat_summary).

## Warning: Removed 327 rows containing non-finite values (stat_summary).

## Warning: Removed 297 rows containing non-finite values (stat_summary).

## Warning: Removed 297 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 6 rows containing non-finite values (stat_summary).

## Warning: Removed 6 rows containing non-finite values (stat_summary).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

## Warning: Removed 8 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 16 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

## Warning: Removed 27 rows containing non-finite values (stat_summary).

These are barplots by GFA of the SMA by quartiles.

# This is based on tabledata, which was generated by the standard deviation variable
# of the GFAs
# Pick the color based on: https://htmlcolorcodes.com/color-picker/

tabledata <- currdata[,c(myqGFA,listvars)]

tablenames <- names(tabledata)
mysmas <- tablenames[grep("_y_",tablenames)]
mysmalabels <- c("Watching TV/Movies per week [h]","Watching Videos per week [h]","Playing Games per week [h]","Texting per week [h]","Social Networking per week [h]","Chatting per week [h]")

for (i in 1:8){
  p <- list()
  for (j in 1: length(mysmas)){
  myplotGFA <- paste0("qGFA",i)
completetabledata <- tabledata[complete.cases(tabledata[,c(myplotGFA)]),]

p[[j]] <-ggplot(data=completetabledata, aes_string(x=myplotGFA, y=mysmas[j])) +
  geom_bar(stat = "summary", fun.y = "mean",fill="#FF3396") + coord_flip() + 
  # ggtitle("Media Activity and GFA") +
  ylim(0, 22) +
  stat_summary(aes(label=round(..y..,2)), fun.y=mean, geom="text", size=3,
             color = "black",hjust=-0.1) +
  xlab(paste0("GFA",i)) + ylab(mysmalabels[j])  + theme_minimal()
# print(p[[i]])
  }
  do.call("grid.arrange", c(p,list(ncol=2, top=textGrob(paste0("GFA ",i)))))
# do.call(grid.arrange,c(p,main="text"))
}

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

for (j in 1: length(mysmas)){
p <- list()
for (i in 1:4){
  myplotGFA <- paste0("qGFA",i)
completetabledata <- tabledata[complete.cases(tabledata[,c(myplotGFA)]),]

p[[i]] <-ggplot(data=completetabledata, aes_string(x=myplotGFA, y=mysmas[j])) +
  geom_bar(stat = "summary", fun.y = "mean",fill="skyblue") + coord_flip() + 
  # ggtitle("Media Activity and GFA") +
  ylim(0, 22) +
  stat_summary(aes(label=round(..y..,2)), fun.y=mean, geom="text", size=3,
             color = "black",hjust=-0.1) +
  xlab(paste0("GFA",i)) + ylab(mysmalabels[j])  + theme_minimal()
# print(p[[i]])
}
do.call("grid.arrange", c(p,list(ncol=2, top=textGrob(paste0("GFA Quartiles and ",mysmalabels[j])))))
# do.call(grid.arrange,p)
}

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 403 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 343 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 359 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 47 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 22 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).

## Warning: Removed 37 rows containing non-finite values (stat_summary).